AD 


Award  Number:  W81XWH-06-1-0003 


TITLE:  Intra-Operative  Dosimetry  in  Prostate  Brachytherapy 

PRINCIPAL  INVESTIGATOR:  AmeetJain 

CONTRACTING  ORGANIZATION:  Johns  Hopkins  University 

Baltimore,  MD  21218-2686 

REPORT  DATE:  November  2007 

TYPE  OF  REPORT:  Annual  Summary 

PREPARED  FOR:  U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


DISTRIBUTION  STATEMENT:  Approved  for  Public  Release; 

Distribution  Unlimited 


The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and 
should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision 
unless  so  designated  by  other  documentation. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

01-11  -2007  Annual  Summary 

3.  DATES  COVERED 

15  Oct  2005 -14  Oct  2007 

4.  TITLE  AND  SUBTITLE 

Intra-Operative  Dosimetry  in  Prostate  Brachytherapy 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

W81XWH-06-1-0003 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Ameet  Jain 

E-Mail:  jain@cs.jhu.edu 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

Johns  Hopkins  University 
Baltimore,  MD  21218-2686 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 


13.  SUPPLEMENTARY  NOTES 

Original  contains  colored  plates:  ALL  DTIC  reproductions  will  be  in  black  and  white. 


14.  ABSTRACT 

Favorable  outcome  in  prostate  brachytherapy  critically  depends  on  the  accurate  placement  of  radioactive  sources  in  their 
planned  locations.  Unfortunately,  there  is  variety  of  mechanical  factors  that  cause  the  seeds  to  divert  from  their  planned 
locations.  While  this  problem  has  been  known  to  brachytherapists,  current  technology  does  not  allow  for  reliable  localization  of 
the  implanted  sources,  thereby  prohibiting  the  prediction  and  modification  of  seed  distribution  intra-operatively.  The  Research 
Objective  of  the  proposal  is  to  develop  and  evaluate  ex-vivo  a  method  for  intra-operative  localization  of  the  implanted  seeds  in 
relation  to  the  prostate,  to  facilitate  in-situ  dosimetric  optimization  and  exit  dosimetry.  In  particular,  we  will:  [1]  Registration  of 
Ultrasound  to  Fluoroscopy  (RUF):  Develop  methods  for  reconstruction  of  seed  implants  from  X-ray  fluoroscopy  and  spatially 
register  them  to  the  prostate  anatomy  identified  in  TRUS  [2]  System  Integration:  Integrate  the  above  methods  in  a  software 
package  and  link  it  with  the  FDA-approved  CMS  Interplant®  prostate  brachytherapy  system  to  enable  in-situ  dosimetry 
calculation  [3]  Experimental  Validation:  Evaluate  the  performance  of  the  RUF  system  on  phantoms  and  pre-recorded  patient 
data. 


15.  SUBJECT  TERMS 

Prostate  Brachytherapy,  X-ray  reconstruction,  C-arm,  TRUS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

USAMRMC 

a.  REPORT 

u 

b.  ABSTRACT 

u 

c.  THIS  PAGE 

u 

uu 

69 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Table  of  Contents 


Introduction . 4 

Body . 4 

Key  Research  Accomplishments . 7 

Reportable  Outcomes . 7 

Conclusions . 7 

References . 8 


Appendices 


9 


Intra-operative  Dosimetry  in  Prostate  Brachytherapy,  PC050170 


Ameet  Jain 


Progress  Report  Summary 

Project  Year  2  (2006  -2007) 


A  INTRODUCTION 

For  several  decades,  the  definitive  treatment  for  low  and  medium  risk  prostate  cancer  was  radical  prostatectomy 
or  external  beam  radiation  therapy,  but  low  dose  rate  permanent  seed  brachytherapy  (shortly  brachytherapy) 
today  can  achieve  equivalent  outcomes.  Brachytherapy,  if  accurately  executed,  can  achieve  a  sharp  demarcation 
between  the  treated  volume  and  healthy  structures,  and  thereby  achieve  superior  tumor  control  with  reduced 
morbidity.  In  contemporary  practice,  however,  faulty  needle  and  source  placement  often  cause  insufficient  dose 
to  the  cancer  and/or  inadvertent  radiation  of  the  rectum,  urethra,  and  bladder.  Another  fallacy  of  the  current 
implant  techniques  is  that  reliable  and  accurate  exit  dosimetry  is  not  possible.  The  contribution  of  the  proposed 
research  will  be  making  C-arm  fluoroscopy  available  for  safe,  simple,  and  robust  intra-operative  localization  of 
brachytherapy  sources  relative  to  the  prostate.  We  will  develop  a  method  for  the  registration  of  ultrasound  to 
fluoroscopy  (RUF),  to  fuse  TRUS  (Trans-rectal  ultrasound  can  view  the  prostate  but  not  the  seeds)  with  C-arm 
fluoroscopy  (which  is  capable  of  viewing  the  seeds  but  not  the  prostate).  This  feature  will  allow  for  dosimetric 
optimization  of  the  prostate  brachytherapy  implants  and  exit  dosimetry  before  the  patient  is  released  from  the 
operating  room;  thereby  enabling  significant  improvement  on  current  clinical  practice.  A  further  promise  is  that 
fluoroscopy-based  exit  dosimetry  may  obviate  CT-based  post-implant  dosimetry. 


B  BODY 

B.l  Brief  System  Concept 

The  system  concept  for  registration  of  ultrasound  to 
fluoroscopy  (RUF)  is  summarized  in  Figure  1.  The 
fluoroscope  is  calibrated  and  corrected  from  image 
distortion  pre-operatively.  The  implant  procedure  starts  as 
usual:  TRUS  is  used  to  guide  each  individual  needle  and  a 
C-arm  placed  over  the  patient’s  abdomen.  The  C-arm  is 
tracked  with  an  X-ray  fiducial  system  called  FTRAC  that  is 
composed  of  optimally  selected  polynomial  space  curves. 

The  fiducials  are  mounted  rigidly  to  the  TRUS  frame  in  the 
field  of  view  in  a  known  calibrated  pose  relative  to  the  TRUS,  thereby  providing  spatial  registration  between  the 
C-arm  and  TRUS.  Upon  implanting  a  batch  of  needles  (typically  a  row  of  needles),  we  collect  a  set  of  TRUS 
and  C-arm  images.  The  locations  of  the  implanted  seeds  are  recovered  from  the  C-arm  fluoroscopy  images  with 
the  use  of  a  network  flow  based  method  called  MARSHAL  that  provides  seed  segmentation,  matching,  and 
reconstruction  method.  Then  the  cloud  of  seeds  is  superimposed  over  the  spatially  co-registered  TRUS  images. 
The  3D  dose  distribution  is  rapidly  calculated  from  the  union  of  the  already  and  yet  to  be  implanted  seeds.  The 
dose  distribution  is  analyzed  with  tools  currently  available  in  the  brachytherapy  system  used.  Then  the  implant 
plan  can  be  optimized  to  account  for  discrepancies  from  the  ideal  dose  distribution.  The  procedure  continues 
with  the  next  batch  of  needles  in  the  cycle  described  above.  After  the  last  needle,  a  complete  dosimetry  check  is 
performed,  which  provides  a  final  opportunity  to  patch  up  any  cold  spots  with  additional  seeds. 

B.2  Proposed  Statement  of  Work 

We  proposed  to  develop  a  method  for  the  registration  of  ultrasound  to  fluoroscopy  (RUF)  to  allow  for  intra¬ 
operative  dosimetry  in  prostate  brachytherapy  and  prototype  mathematical  algorithms  (Aim-1),  integrate  them 
with  an  existing  FDA  approved  prostate  brachytherapy  system  that  provides  dosimetry  analysis  (Aim-2),  and 
evaluate  the  system  experimentally  on  phantoms  and  pre-recorded  patient  data  (Aim-3).  Algorithmic  design 
(Aim-1)  and  experimental  evaluation  (Aim-3),  will  progress  hand  in  hand.  System  integration  (Aim-2)  will  be 
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performed  immediately  when  a  workable  subset  of  RUF  package  becomes  available  from  Aim- 1  and  again 
revisited  towards  at  the  end  of  the  project.  Therefore,  the  timeline  will  be  somewhat  non-linear.  The  detailed 
statement  of  work  was  as  follows: 

Aim-1:  Registration  of  Ultrasound  to  Fluoroscopy  (RUF):  Develop  a  methods  for  reconstruction  of  seed 
implants  from  X-ray  fluoroscopy  and  spatially  registering  them  to  the  prostate  anatomy  identified  in 
TRUS 

Aim-2:  System  Integration:  Integrate  the  above  methods  in  a  software  package  and  link  it  with  the  FDA- 
approved  CMS  Interplant®  prostate  brachytherapy  system  to  enable  dosimetry  calculation 
Aim-3:  Experimental  Validation:  Evaluate  the  performance  of  the  RUF  system  on  phantoms  and  pre¬ 
recorded  patient  data.  (Neither  of  which  require  an  IRB  approval) 

B.3  Progress  Report  for  Second  Year 

In  Aiml:  the  three  main  aspects  were  (a)  Registration  of  X-ray  to  ultrasound;  (b)  C-arm  tracking;  (c)  Seed 
reconstruction.;  and  (d)  C-arm  calibration. 

These  issues  were  satisfactorily  addressed  in  the  year-1  report. 


In  Aim2:  we  had  reported  preliminary  work  and  proposed  approach  towards  the  system  integration  in  the  year  1 
annual  report.  We  made  further  progress  and  completed  this  integration  process  in  the  second  year. 


We  employ  a  regular  clinical  brachytherapy  setup,  without  alteration,  including  a  treatment  planning 
workstation  and  stabilizer/stepper  (Interplant®,  CMS,  St  Louis),  TRUS  (B&K  Medical),  C-arm  (GE  OEC,  Salt 
Lake  City,  UT).  The  C-arm  is  interfaced  with  a  standalone  laptop  computer  through  an  NTSC  video  line  and 
frame  grabber. 


Workflow:  The  clinical  workflow  is  identical  to  the  standard  procedure  until  the  clinician  decides  to  run  a 
reconstruction  and  optimization,  typically  at  each  half  of  the  procedure.  A  set  of  C-arm  images  are  collected 
with  as  wide  separation  as  possible  (15°)  and  transferred  to  the  laptop.  After  processing  the  fluoro  images,  the 
seeds  are  reconstructed  in  TRUS  space  and  exported  to  the  Interplant®.  The  physician  uses  standard 
Interplant®  features  to  analyze  the  dose  and  modify  the  remainder  of  the  implant  plan,  in  order  to  optimize  it. 
The  procedure  concludes  when  the  exit  dosimetry  shows  no  cold  spots. 


Dosimetry  Analysis  and  Implant  Optimization:  The  seed  locations  are  exported  in  template  coordinates  to  the 
Interplant®  system.  A  software  patch  added  to  the  Interplant®  then  removes  the  already  implanted  seeds  from 
the  complete  implant  plan,  thereby  producing  a  “residual  implant  plan”.  The  total  dose  is  calculated  by 
combining  the  seeds  in  the  residual  implant  plan  and  the  seeds  already  implanted.  The  physician  uses  standard 
Interplant®  tools  (isodose  coverage,  DVH,  etc.)  for  dose  analysis.  If  it  is  necessary  to  change  the  total  dose 
distribution,  then  the  residual  plan  is  modified  with  standard  Interplant®  dose  planning  functions.  Note  that  as 
the  procedure  progresses,  the  residual  plan  reduces  and  there  is  less  and  less  degree  of  freedom  to  optimize  the 
overall  dose  distribution.  In  fact,  at  some  point  it  becomes  impossible  to  reduce  dose  in  the  hot  spots,  because 
once  seeds  are  delivered,  they  cannot  be  taken  back.  Fortunately,  we  always  can  fix  cold  spots  by  adding  more 
seeds  to  the  implant  plan. 


In  Aim3:  we  proposed  to  evaluate  the  workspace  constraints  and  the  overall  performance  of  the  system.  In 
addition  to  the  results  presented  at  the  end  of  year  1,  further  extensive  validation  was  conducted,  including 
detailed  phantom  experiments  for  engineering  validation  and  some  clinical  validation. 
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Soft  Training  Phantoms:  We  fully  seeded  three  standard  prostate  brachytherapy  phantoms  with  realistic  implant 
plans  (45,  49,  87  seeds).  Seed  locations  reconstructed  from  fluoro  using  a 
realistic  15°  net  C-arm  motion  were  compared  to  corresponding  locations 
segmented  manually  in  CT  (1mm  slice  thickness).  The  average  seed 
reconstruction  error  wrt  the  FTRAC  fiducial  was  below  1.5mm  (0.2mm 
relative  accuracy  -  there  was  a  translation  offset).  These  results  made  us 
confident  that  the  system  would  perform  reasonably  for  actual  humans.  To 
measure  the  full  system  error,  5  needles  (tips)  were  inserted,  reconstructed 
and  exported  into  a  standard  prostate  brachytherapy  phantom.  Manual 
segmentation  of  the  needles  in  TRUS  (sagittal  for  depth  and  transverse  for 
X-Y)  provided  ground  truth  locations.  The  average  residual  error  for  the  5 

needle  tips  was  3.4  mm.  There  is  about  a  3mm  constant  translation  Figure  2:  Prostate  phantom  with  the  C- 
offset/bias,  due  to  cascading  effects  in  the  Interplant®  pre-operative 


arm,  FTRAC,  TRUS  and  the  Interplant® 


registration,  the  reconstruction  process  &  US  beam  thickness.  After  correcting  for  these  biases,  the  mean 
residual  registration  error  dropped  to  0.82mm  (STD=0.17mm). 


Clinical  Data:  From  a  small  cohort  of  6  patients  so  far  we  have  treated  4.  An  actual  OR  picture  in  Figure  2 
shows  the  setup  with  the  FTRAC  fluoroscopy 
tracking  fiducial  mounted  over  the  template  and 
TRUS  stepper. 

Implant  reconstruction  and  dose  analysis  was 
performed  halfway  of  the  implant  and  also  in 
the  end  (an  additional  one  was  used  for  a 
research  comparison  to  post-op  CT).  In  all  the 
patients,  the  final  dosimetry  detected  one  or 
more  cold  spots  (Figure  5).  In  the  fourth  patient, 

9  additional  seeds  (initial  67)  were  inserted  to 
bring  up  VI 00  dose  to  the  desired  level. 

Apparently,  the  physician  grew  quickly  to  trust 


Figure  3:  Fluoroscope  tracking  (FTRAC)  fiducial:  Mounting  over  the 


abdomen  in  CAD  model  (left)  and  actual  clinical  setup  (middle),  the 
connector  between  FTRAC  and  template  (right). 


ip 


.ip  I _ l _ I _ ; _ l _ 1 

lu  p  £  M  IS  3}  3S  3P  36 

Figure  4:  Intra-operative  visualization  of  the  edema  is  possible  with 
the  system.  The  planned  (green)  versus  the  ‘ reconstructed ’  (red)  seed 
positions  as  seen  in  the  X-Y  plane  (template  view).  The  image 
acquisition  is  at  end  of  procedure.  A  trend  of  outward  radiation  of  the 
seeds  from  their  initial  seed  position  is  observed.  This  is  due  to  edema 
occurring  during  the  procedure. 


the  system  in  catching  cold  spots,  and  instead 
became  very  cautions  not  to  create  hot  spots. 
Importantly,  all  patients  were  released  from  the 
operating  room  with  satisfactory  V  I 00/V90  dose 
and  favorable  DVH  for  all  relevant  structures.  The 
use  of  the  system  added  about  45  minutes  to  the 
OR  time,  with  each  set  of  three  implant 
reconstructions  taking  about  15  minutes, 
including  image  capture,  processing,  3D 
reconstruction,  and  dosimetry  evaluation.  We 
anticipate  needing  only  a  single  reconstruction  in 
the  post-research  phase.  Intra-operative 
visualization  of  edema  was  also  possible  with  the 
system,  which  was  found  to  be  significant  in  2 
patients  of  the  total  4  (Figure  4).  The  seeds 
showed  a  clear  tendency  of  outward  migration 
from  their  drop  positions  as  the  procedure 
progressed.  Although  comparative  quantitative 
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conclusions  cannot  be  made  based  on  the  limited  available  data,  comparison  of  the  intra-operative  dosimetry 
with  the  post-implant  CT  reveal  close  correlation. 


Figure  5:  The  system  is  able  to  detect  cold  spots.  The 
VI 00  contours  (blue)  as  assumed  by  the  planning  system 
(top)  and  as  computed  by  the  proposed  system  (bottom).  In 
this  patient,  3  cold  spots  were  discovered  in  exit  dosimetry. 
The  prostate  boundary  is  delineated  by  dark  red  contours 
and  the  green  squares  mark  current  seed  coordinates.  Note 
that  the  system  also  detected  4  seeds  that  have  drifted  out 
of  slice,  owing  largely  to  the  swelling  effects  of  edema. 


C  KEY  RESEARCH  ACCOMPLISHMENTS 

1 .  Ability  to  intra-operative ly  compute  3D  seed  locations  in  prostate  brachytherapy. 

2.  Minimal  alterations  to  the  current  prostate  brachytherapy  clinical  protocol. 

3.  A  complete  integrated  clinical  research  system  ready  for  multi-center  trials. 

4.  Extensive  validation  -  both  on  phantom  data  as  promised,  and  additionally  also  on  real  clinical  data. 

D  KEY  TRAINING  ACCOMPLISHMENTS 

1 .  Close  interaction  with  clinicians,  medical  physicists,  and  industry  advisors. 

2.  Various  international  training  courses,  workshops  &  conferences  in  medical  imaging  processing. 

3.  Ability  to  do  independent  research,  including  mentoring  undergraduate  students  for  their  research. 

4.  Successful  completion  of  the  PhD  requirements. 

E  REPORTABLE  OUTCOMES 

1.  Publications  in  leading  conferences  and  journals,  as  listed  in  the  publications  section  (including  three 
flagship  conference  publications  in  MICCAI). 

2.  A  poster  award  at  SPIE  Medical  Imaging  2007  for  work  on  C-arm  calibration. 

3.  A  prototype  for  subsequent  commercialization  and  integration  into  the  Interplant  software. 

F  CONCLUSIONS 

In  conclusion,  the  system  was  exquisitely  accurate  in  technical  (phantom)  trials  and  has  shown  usefulness  and 
great  potential  in  a  limited  trial.  Further  development  is  underway  to  reduce  added  time  and  manual  labor  in  the 
OR.  We  have  addressed  the  issue  of  intra-operative  reconstruction  of  brachytherapy  seeds,  with  minimal 
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alteration  to  the  current  clinical  protocol  or  any  significant  increase  in  cost.  A  prototype  system  for  clinical  trails 
and  potential  commercialization  has  also  been  tested  &  implemented. 


So  what:  The  success  of  brachytherapy  chiefly  depends  on  our  ability  to  intra-operatively  cover  the  prostate 
with  sufficient  radiation  while  still  avoiding  excessive  radiation  to  surrounding  organs.  Currently,  such  level  of 
precision  is  not  always  achievable  even  by  the  most  experienced  physicians.  Thus  many  implants  fail  or  cause 
severe  side  effects  owing  to  faulty  seed  placement,  a  problem  what  still  cannot  be  corrected  in  the  operating 
room.  Our  results  indicate  the  feasibility  of  a  system  that  could  achieve  intra-operative  localization  of  the 
implanted  seeds  in  relation  to  the  prostate,  to  allow  for  in-situ  dosimetric  optimization  and  exit  dosimetry.  This 
ability  to  perform  intra-operative  dosimetry  may  change  the  standard  of  care  in  brachytherapy  by  allowing  the 
physician  to  achieve  technically  excellent  brachytherapy  implants,  resulting  in  improved  disease  control  and 
quality  of  life  for  a  large  and  steadily  growing  group  of  patients. 

G  PUBLICATIONS  (List  of  publications  in  year  2  arising  from  the  grant,  sorted  by  publication  date) 
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Intra-operative  Guidance  in  Prostate  Brachytherapy  Using  an  non-isocentric  C-arm.  Tenth  International 
Conference  on  Medical  Image  Computing  and  Computer- Assisted  Intervention  (MICCAI),  Oct  2007 
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Abstract.  Intra-operative  guidance  in  Transrectal  Ultrasound  (TRUS) 
guided  prostate  brachytherapy  requires  localization  of  inserted  radioac¬ 
tive  seeds  relative  to  the  prostate.  Seeds  were  reconstructed  using  a 
typical  C-arm,  and  exported  to  a  commercial  brachytherapy  system  for 
dosimetry  analysis.  Technical  obstacles  for  3D  reconstruction  on  a  non- 
isocentric  C-arm  included  pose-dependent  C-arm  calibration;  distortion 
correction;  pose  estimation  of  C-arm  images;  seed  reconstruction;  and 
C-arm  to  TRUS  registration.  In  precision-machined  hard  phantoms  with 
40-100  seeds,  we  correctly  reconstructed  99.8%  seeds  with  a  mean  3D  ac¬ 
curacy  of  0.68  mm.  In  soft  tissue  phantoms  with  45-87  seeds  and  clinically 
realistic  15°  C-arm  motion,  we  correctly  reconstructed  100%  seeds  with 
an  accuracy  of  1.3  mm.  The  reconstructed  3D  seed  positions  were  then 
registered  to  the  prostate  segmented  from  TRUS.  In  a  Phase- 1  clinical 
trial,  so  far  on  4  patients  with  66-84  seeds,  we  achieved  intra-operative 
monitoring  of  seed  distribution  and  dosimetry.  We  optimized  the  100% 
prescribed  iso-dose  contour  by  inserting  an  average  of  3.75  additional 
seeds,  making  intra-operative  dosimetry  possible  on  a  typical  C-arm,  at 
negligible  additional  cost  to  the  existing  clinical  installation. 


1  Introduction 

With  an  approximate  annual  incidence  of  220,000  new  cases  and  33,000  deaths 
(United  States)  prostate  cancer  continues  to  be  the  most  common  cancer  in  men. 
Transrectal  Ultrasound  (TRUS)  guided  permanent  low-dose-rate  brachytherapy 
(insertion  of  radioactive  seeds  into  the  prostate)  has  emerged  as  a  common  & 
effective  treatment  modality  for  early  stage  low  risk  prostate  cancer,  with  an 
expected  50,000  surgeries  every  year.  The  success  of  brachytherapy  (i.e.  max¬ 
imizing  its  efficacy  while  minimizing  its  co-morbidity)  chiefly  depends  on  our 
ability  to  tailor  the  therapeutic  dose  to  the  patient’s  individual  anatomy.  The 
main  limitation  in  contemporary  brachytherapy  is  intra-operative  tissue  expan¬ 
sion  (edema),  causing  incorrect  seed  placement,  which  may  potentially  lead  to 
insufficient  dose  to  the  cancer  and/or  excessive  radiation  to  the  rectum,  ure¬ 
thra,  or  bladder.  The  former  might  permit  the  cancer  to  relapse,  while  the  latter 
causes  adverse  side  effects  like  rectal  ulceration.  According  to  a  comprehensive 
review  by  the  American  Brachytherapy  Society  [1],  the  pre-planned  technique 
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used  for  permanent  prostate  brachytherapy  has  limitations  that  may  be  overcome 
by  intra-operative  planning. 

Prostate  brachytherapy  is  almost  exclusively  performed  under  TRUS  guid¬ 
ance.  Various  researchers  have  tried  to  segment  the  seeds  from  TRUS  images 
by  linking  seeds  with  spacers,  using  X-rays  to  initialize  segmentation,  using 
vibro-acoustography  or  transurethral  ultrasound  as  a  new  imaging  modality, 
or  segmenting  them  directly  in  TRUS  images  by  using  corrugated  seeds  that 
are  better  visible  than  conventional  ones  [2] .  But  even  when  meticulously  hand- 
segmented,  up  to  25%  of  the  seeds  may  remain  hidden  in  ultrasound.  C-arms  are 
also  ubiquitous,  though  used  only  for  gross  visual  assessment  of  the  implanted 
seed  positions  (approximately  60%  of  the  practitioners  using  it  in  the  operating 
room  [3]).  In  spite  of  significant  efforts  that  have  been  made  towards  computa¬ 
tional  fluoroscopic  guidance  in  general  surgery  [4],  C-arms  cannot  yet  be  used 
for  intra-operative  brachytherapy  guidance  due  to  a  plethora  of  technical  lim¬ 
itations.  While  several  groups  have  published  protocols  and  clinical  outcomes 
favorably  supporting  C-arm  fluoroscopy  for  intra-operative  dosimetric  analysis 
[5-7],  this  technique  is  yet  to  become  a  standard  of  care  across  hospitals.  In  this 
paper  we  report  a  system  to  reconstruct  3D  seed  positions  (visible  in  X-ray)  and 
spatially  register  them  to  the  prostate  (visible  in  TRUS).  Our  primary  con¬ 
tribution  is  our  ability  to  use  any  typical  non-isocentric  uncalibrated  C-arm 
present  in  most  hospitals,  in  comparison  to  the  use  of  calibrated  isocentric  ma¬ 
chines  [5,6]  or  an  approximate  reconstruction  [7],  as  reported  in  the  literature. 


2  Methods  and  Materials 

The  system  is  designed  to  easily  integrate  with  commercial  brachytherapy  instal¬ 
lations.  We  employ  a  regular  clinical  brachytherapy  setup,  without  alteration, 
including  a  treatment  planning  workstation  &  stabilizer/stepper  (Interplant®, 
CMS,  St  Louis),  TRUS  (B&K  Medical  Pro  Focus)  and  a  C-arm  (GE  OEC 
98/9900).  The  C-arm  is  interfaced  with  a  laptop  through  an  NTSC  video  line 
and  frame  grabber,  making  the  image  capture  independent  of  the  C-arm  model. 
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Fig.  1.  Overview  of  the  proposed  solution.  The  FTRAC  fiducial  tracks  C-arms,  and 
also  registers  TRUS  to  C-arm  images,  making  quantitative  brachytherapy  possible. 


Workflow:  The  clinical  workflow  (Fig.  1)  is  identical  to  the  standard  procedure 
until  the  clinician  decides  to  run  a  reconstruction  and  optimization.  A  set  of  C- 
arm  images  are  collected  with  a  separation  as  wide  as  clinically  possible  (10  —  15° 
around  AP-axis)  and  synchronously  transferred  to  the  laptop.  After  processing 
the  images,  the  seeds  are  reconstructed  and  exported  to  the  Interplant®  system. 
The  physician  uses  standard  Interplant®  tools  to  analyze,  optimize  and  modify 
the  remainder  of  the  plan.  The  procedure  concludes  when  the  exit  dosimetry 
shows  no  cold  spots  (under-radiated  locations). 

Numerous  technical  obstacles  have  to  be  overcome  to  realize  C-arm  based  intra¬ 
operative  dosimetry:  (a)  C-arm  calibration;  (b)  image  distortion  correction;  (c) 
pose  estimation  of  C-arm  images;  (d)  seed  reconstruction;  (e)  registration  of  C- 
arm  to  TRUS;  (f)  dosimetry  analysis;  and  finally  (g)  implant  optimization.  We 
have  developed  a  system  that  overcomes  these  limitations  in  providing  quanti¬ 
tative  intra-operative  dosimetry.  In  what  follows,  we  will  describe  briefly  each 
component  of  the  system,  skipping  the  mathematical  framework  for  lack  of  space. 

C-arm  Source  Calibration  and  Image  Distortion:  Since  both  C-arm  cali¬ 
bration  and  distortion  are  pose-dependent,  contemporary  fluoroscopy  calibrates/ 
distortion-corrects  at  each  imaging  pose  using  a  cumbersome  calibration- fixture, 
which  is  a  significant  liability.  Our  approach  is  a 
complete  departure.  Using  a  mathematical  &  ex¬ 
perimental  framework,  we  demonstrated  that  cali¬ 
bration  is  not  critical  for  prostate  seed  reconstruc¬ 
tion.  Just  an  approximate  pre-operative  calibra¬ 
tion  suffices  [8].  The  central  intuition  is  that  ob¬ 
ject  reconstruction  using  a  mis-calibrated  C-arm 
changes  only  the  absolute  positions  of  the  objects, 
but  not  their  relative  ones  (Fig.  2).  Additionally, 
statistical  analysis  of  the  distortion  in  a  15°  lim¬ 
ited  workspace  around  the  AP-axis  revealed  that 
just  a  pre-operative  correction  can  reduce  the  aver¬ 
age  distortion  from  3.31  mm  to  0.51  mm,  sufficient 
for  accurate  3D  reconstruction.  The  numbers  are 
expected  to  be  similar  for  other  C-arms  too. 

Pose  Estimation:  The  most  critical  component  for  3D  reconstruction  is  C-arm 
pose  estimation.  C-arms  available  in  most  hospitals  do  not  have  encoded  rota¬ 
tional  joints,  making  the  amount  of  C-arm  motion  unavailable.  C-arm  tracking 
using  auxiliary  trackers  is  expensive,  inaccurate  in  the  presence  of  metal  (EM 
tracking)  or  intrudes  in  the  operating  room  (optical  tracking).  There  has  been 
some  work  on  fiducial  based  tracking,  wherein  a  fiducial  (usually  large  for  accu¬ 
racy)  is  introduced  in  the  X-ray  FOV  and  its  projection  in  the  image  encodes  the 
6  DOF  pose  of  the  C-arm.  We  proposed  a  new  fluoroscope  tracking  (FTRAC) 
[9]  fiducial  design  that  uses  an  ellipse  (key  contribution),  allowing  for  a  small 
(3x3x5cm)  yet  accurate  fiducial.  In  particular,  the  small  size  makes  it  easier  to 


Fig.  2.  Mis-calibration 
conserves  relative  recon¬ 
struction  between  objects 
A  and  B  (eg.  seeds). 


be  always  in  the  FOV  &  to  be  robust  to  image  dis¬ 
tortion.  Extensive  phantom  experiments  indicated 
a  mean  tracking  accuracy  on  distorted  C-arms  of 
0.56  mm  in  translation  and  0.25°  in  rotation,  an  ac¬ 
curacy  comparable  to  expensive  external  trackers. 

Seed  Segmentation:  We  developed  an  automated 
seed  segmentation  algorithm  that  employs  the  mor¬ 
phological  top-hat  transform  to  perform  the  basic 
seed  segmentation,  followed  by  thresholding,  region 
labeling,  and  finally  a  two-phase  classification  to  seg¬ 
ment  both  single  seeds  &  clusters.  The  result  of  the 
segmentation  is  verified  on  the  screen  to  allow  for  a 
manual  bypass  by  the  surgeon. 


Fig.  3.  The  FTRAC 
fiducial  mounted  over 
the  seed-insertion  needle 
template  using  a  me¬ 
chanical  connector.  An 
X-ray  image  of  the  same. 


Seed  Correspondence  &  Reconstruction:  The  3D  coordinates  of  the  im¬ 
planted  seeds  can  now  be  triangulated  by  resolving  the  correspondence  of  seeds 
in  the  multiple  X-ray  images.  We  formalized  seed  correspondence  to  a  network- 
flow-based  combinatorial  optimization,  wherein  the  desired  solution  is  the  flow 
with  minimum  cost.  Using  this  abstraction,  we  proposed  an  algorithm  (MAR¬ 
SHAL  [10])  that  runs  in  cubic-time  using  any  number  of  images.  In  comparison, 
previous  solutions  have  predominantly  been  heuristic  explorations  of  the  large 
search  space  (lO300).  In  addition,  the  framework  also  robustly  resolves  all  the 
seeds  that  are  hidden  in  the  images  (typically  4-7%  due  to  the  high  density). 
MARSHAL  typically  reconstructs  99.8%  of  the  seeds  and  runs  in  under  5s  in 
MATLAB  (a  95%  minimum-detection-rate  is  usually  deemed  sufficient  [11]). 

Registration  of  C-arm  to  TRUS:  The  FTRAC  is 

attached  to  the  needle-insertion  template  by  a  pre¬ 
cisely  machined  mechanical  connector  (Fig.  4)  in  a 
known  relative  way  (pre-calibration).  The  template 
has  already  been  calibrated  to  TRUS  as  per  the  usual 
clinical  protocol.  Thus  a  simple  application  of  the 
various  known  frame  transformations,  registers  the 
3D  seeds  (FTRAC)  to  the  prostate  (TRUS). 

System  Implementation ,  Dosimetry  Analysis 
and  Implant  Optimization:  We  have  integrated 
all  the  above  functions  into  a  MATLAB  program  with  a  GUI.  The  package  runs 
on  a  laptop  that  sends  reconstructed  seed  positions  (in  template  coordinates) 
to  the  Interplant®  system.  In  order  to  not  require  a  new  FDA  approval,  we 
maintain  the  integrity  of  the  FDA-approved  Interplant®  by  not  modifying  the 
commercial  software,  and  instead  use  a  text  file  to  export  the  3D  seed  locations. 
The  physician  uses  standard  Interplant®  tools  (isodose  coverage,  etc.)  for  dose 
analysis,  and  if  needed,  modifies  the  residual  plan  to  avoid  hot  spots  or  fill  in 
cold  spots.  This  process  can  be  repeated  multiple  times  during  the  surgery. 


Fig.  4.  FTRAC  &  tem¬ 
plate  pre-calibration  us¬ 
ing  a  rigid  mount. 


3  Phantom  Experiments  and  Results 


We  have  extensively  tested  the  system  and  its  components  in  various  phantoms 
and  in  an  ongoing  Phase-1  clinical  trial.  To  do  so,  we  introduce  the  terms  ab¬ 
solute  and  relative  reconstruction  errors.  Using  X-ray  images,  the  seeds  are 
reconstructed  with  respect  to  (w.r.t.)  the  FTRAC  frame.  In  experiments  where 
the  ground  truth  location  of  the  seeds  w.r.t.  the  FTRAC  is  known,  the  com¬ 
parative  analysis  is  called  absolute  accuracy.  Sometimes  (eg.  in  patients),  the 
true  seed  locations  w.r.t.  the  FTRAC  are  not  available  and  the  reconstruction 
can  only  be  compared  to  the  seeds  extracted  from  post-op  data  (using  a  rigid 
point-cloud  transform),  in  which  case  the  evaluation  is  called  relative  accuracy. 

Solid  Seed  Phantom:  An  acetol  (Delrin)  phantom  consisting  of  ten  slabs 
(5mm  each)  was  fabricated  (Fig.  5  (a)).  This  phantom  provides  a  multitude  of 
implants  with  sub- mm  ground  truth  accuracy.  The  fiducial  was  rigidly  attached 
to  the  phantom  in  a  known  way,  establishing  the  accurate  ground  truth  3D 
location  of  each  seed.  Realistic  prostate  implants  (1.56  seeds/cc,  40-100  seeds) 
were  imaged  within  a  30°  cone  around  the  AP-axis.  The  true  correspondence 
was  manually  established  by  using  the  3D  locations,  known  from  the  precise 
fabrication.  Averaged  results  indicate  that  we  correctly  match  98.5%  &  99.8% 
of  the  seeds  using  3  &  4  images  (100  &  75  total  trials)  respectively.  The  mean 
3D  absolute  reconstruction  accuracy  was  0.65  mm  (STD  0.27  mm),  while  the 
relative  accuracy  was  0.35  mm.  Furthermore,  using  4  images  yielded  only  one 
poorly  mis-matched  seed  from  the  75  datasets,  suggesting  the  use  of  4  images 
for  better  clinical  guidance. 


Fig.  5.  (a)  An  image  of  the  solid  seed  phantom  attached  to  the  fiducial  with  a  typical 
X-ray  image  of  the  combination,  (b)  An  annotated  image  of  the  experimental  setup  for 
the  training  phantom  experiment,  (c)  The  clinical  setup  from  the  Phase-I  clinical  trial 


Soft  Training  Phantoms:  We  fully  seeded  three  standard  prostate  brachy ther¬ 
apy  phantoms  (Fig.  5  (b))  with  realistic  implant  plans  (45,  49,  87  seeds).  Seed 
locations  reconstructed  from  fluoro  using  realistic  (maximum  available  clinically) 
image  separation  (about  15°)  were  compared  to  their  corresponding  ground  truth 
locations  segmented  manually  in  CT  (1mm  slice  thickness).  Additionally,  the  45 
&  87-seed  phantoms  were  rigidly  attached  to  the  FTRAC,  providing  the  absolute 
ground  truth  (from  CT).  The  49-seed  phantom  was  used  to  conduct  a  full  scale 


practice-surgery,  in  which  case  the  3D  reconstruction  could  be  compared  only 
to  the  seed  cloud  from  post-op  CT  (without  FTRAC),  providing  just  relative 
accuracy.  Note  that  our  reconstruction  accuracy  ( as  evident  from  the  previous  ex¬ 
periments)  is  better  than  the  CT  resolution.  The  absolute  reconstruction  errors 
for  the  45,  87-seed  phantoms  were  1.64  mm  &  0.95  mm  (STD  0.17  mm),  while 
the  relative  reconstruction  errors  for  the  45,  49,  87-seed  phantoms  were  0.22  mm, 
0.29  mm,  0.20  mm  (STD  0.13  mm).  A  mean  translation  shift  of  1.32  mm  was 
observed  in  the  3D  reconstructions,  predominantly  due  to  the  limited  C-arm 
workspace  (solid-phantom  experiments  with  30°  motion  have  0.65  mm  accu¬ 
racy)  .  It  was  observed  that  the  shift  was  mostly  random  &  not  in  any  particular 
direction.  Nevertheless,  the  accuracy  is  sufficient  for  br  achy  therapy,  especially 
since  a  small  shift  still  detects  the  cold  spots. 

Patients:  A  total  of  11  batches  of  reconstructions  were  carried  out  on  4  patients 
with  2  —  3  batches/patient  &  22  — 84  seeds/batch.  Since  the  seeds  migrate  by  the 
time  a  post-operative  CT  is  taken,  there  is  no  easy  ground  truth  for  real  patients. 
Hence,  for  each  reconstruction,  5  —  6  additional  X-ray  images  were  taken.  The 
reconstructed  3D  seed  locations  were  projected  on  these  additional  images  and 
compared  to  their  segmented  corresponding  2D  locations.  The  results  from  55 
projections  gave  a  2D  mean  error  of  1.59  mm  (STD  0.33  mm,  max  2.44  mm), 
indicating  a  sub- mm  3D  accuracy  (errors  get  magnified  when  projected). 

Registration  Accuracy:  To  measure  the  accuracy  of  the  fiducial  to  template 
registration,  three  batches  of  five  needles  each  were  inserted  randomly  at  random 
depths  into  the  template.  Their  reconstructed  tip  locations  were  then  compared 
to  their  true  measured  locations  (both  in  template  coordinates).  The  limited- 
angle  image-capture  protocol  was  kept  similar  to  that  used  in  the  clinic.  The 
average  absolute  error  (reconstruction  together  with  registration)  was  1.03  mm 
(STD  0.20  mm),  while  the  average  relative  error  was  0.36  mm  (STD  0.31  mm), 
with  an  average  translation  shift  of  0.97  mm. 

System  Accuracy:  To  measure  the  full  system  error,  5  needles  (tips)  were 
inserted  into  a  prostate  brachytherapy  training  phantom,  reconstructed  in  3D 
and  exported  to  the  Interplant®  software.  Manual  segmentation  of  the  needles  in 
TRUS  images  (sagittal  for  depth  and  transverse  for  X-Y)  provided  ground  truth. 
The  mean  absolute  error  for  the  5  needle  tips  was  4  mm  (STD  0.53  mm),  with  a 
translation  shift  of  3.94  mm.  In  comparison,  the  relative  error  for  the  complete 
system  was  only  0.83  mm  (STD  0.18  mm).  The  shift  can  mainly  be  attributed 
to  a  bias  in  the  Template-TRUS  pre-calibration  (~  3  mm)  done  as  part  of  cur¬ 
rent  clinical  practice,  &  in  the  3D  reconstruction  (^  1  mm).  Nevertheless,  we 
removed  this  shift  in  the  clinical  cases  by  applying  a  translation-offset  to  the  re¬ 
constructed  X-ray  seed  coordinates.  This  offset  was  intra-operatively  estimated 
by  comparing  the  centroid  of  the  reconstructed  seeds  with  that  of  the  planned 
seed  locations,  and  by  aligning  the  two  together.  Note  that  the  centroid  is  a 
first-order  statistic  and  robust  to  any  spatially  symmetric  noise/displacement 
model.  Though  a  heuristic,  it  provided  excellent  qualitative  results  according 


to  the  surgeon,  who  read  the  visual  cues  at  the  reconstructed  seed  locations  in 
TRUS  images.  Based  on  the  experiments  so  far  and  the  surgeon’s  feedback,  the 
overall  accuracy  of  the  system  is  expected  to  be  1  —  2  mm  during  clinical  use. 


(a)  (b) 

Fig.  6.  (a)  The  system  is  able  to  detect  cold  spots.  The  100%  iso-dose  contours  (pink) 
as  assumed  by  the  planning  system  (top)  and  as  computed  by  the  proposed  system 
(bottom),  discovering  2  cold  spots.  Red  marks  the  prostate  boundary.  The  green  squares 
delineate  the  seed  coordinates,  detecting  4  seeds  that  had  drifted  out  of  slice,  (b)  The 
system  can  visualize  intra-operative  edema  (mean  4.6  mm,  STD  2.4  mm,  max  12.3  mm). 
The  ’planned’  (red)  versus  the  ’reconstructed’  (blue)  seed  positions  as  seen  in  the 
template  view.  A  trend  of  outward  radiation  from  their  initial  locations  is  observed. 

Phase- 1  Clinical  Trial:  We  have  treated  4  patients  so  far  (Fig.  5  (c)),  out  of  a 
total  of  6  that  will  be  enrolled.  Intra-operative  dosimetry  was  performed  halfway 
during  the  surgery,  at  the  end,  and  after  additional  seed  insertions.  The  current 
protocol  adds  15  minutes  to  the  OR  time  for  each  reconstruction,  including  the 
capture  of  extra  images  (validation),  reconstruction,  and  dosimetry  optimization. 
In  regular  clinical  practice,  we  anticipate  the  need  for  only  a  single  exit-dosimetry 
reconstruction,  increasing  the  OR  time  by  about  10  minutes.  In  all  the  patients 
the  final  dosimetry  detected  cold  spots  (Fig.  6  (a)).  The  clinician  grew  quickly 
to  trust  the  system  in  detecting  cold  spots,  and  instead  minimized  potential  hot 
spots  during  the  surgery.  All  patients  were  released  from  the  operating  room  with 
satisfactory  outcomes.  Intra-operative  visualization  of  edema  (prostate  swelling) 
was  also  possible  (Fig.  6  (b)),  which  was  found  to  be  0.73,  4.64,  4.59,  4.05  mm 
(STD  1.1,  2.2,  2.34,  2.37  mm).  The  seeds  (and  hence  the  prostate)  showed  a  clear 
tendency  for  outward  migration  from  their  drop  positions  (with  maximums  up 
to  15  mm).  Edema  is  the  single  largest  factor  that  makes  the  perfect  delivery 
of  the  pre-planned  dose  nearly  impossible.  In  almost  all  the  patients,  towards 
the  end  of  the  surgery,  it  was  found  that  the  apex  of  the  prostate  (surgeon-end) 
was  under-dosed.  The  medical  team  found  the  intra-operative  visualization  of 
under-dosed  regions  valuable,  inserting  an  additional  1,  2,  3,  9  seeds  to  make  the 
100%  prescribed  iso-dose  contour  cover  the  prostate.  A  further  comparison  of  the 
exit  implant  to  Day-0  CTs  (2  mm  slices)  showed  mean  errors  of  5.43,  6.16,  3.13, 


5.15  mm  (STD  2.46,  2.96,  2.02,  2.71  mm),  indicating  a  further  post-operative 
seed  migration.  Though,  post-operative  seed  migration  is  an  inherent  limitation 
in  br  achy  therapy,  surgeons  usually  accommodate  for  it  by  slightly  over-dosing 
the  patient  (note  that  sub- mm  seed  placement  is  non-critical).  A  study  with  40 
patients  is  currently  being  planned,  to  make  a  statistically  relevant  evaluation 
of  the  medical  benefit  of  the  system  using  clinical  indicators. 

4  Conclusion,  Shortcomings  and  Future  Work 

A  system  for  brachytherapy  seed  reconstruction  has  been  presented,  with  exten¬ 
sive  phantom  and  clinical  trials.  The  absolute  seed  reconstruction  accuracy  from 
phantom  trials  is  1.3  mm  using  15°  C-arm  motion,  sufficient  for  detection  of 
any  cold  spots.  It  shows  usefulness  and  great  potential  from  the  limited  Phase- 1 
patient  trials.  The  system  (a)  requires  no  significant  hardware;  (b)  does  not  alter 
the  current  clinical  workflow;  (c)  can  be  used  with  any  C-arm;  &  (d)  integrates 
easily  with  any  pre-existing  brachytherapy  installation,  making  it  economically 
sustainable  and  scalable.  There  is  some  added  radiation  to  the  patient,  though 
insignificant  when  compared  to  that  from  the  seeds.  Though  not  critical,  pri¬ 
mary  shortcomings  include  (a)  15  minute  additional  OR  time;  (b)  supervision 
during  segmentation;  &  (c)  a  small  translation  bias.  Furthermore,  a  TRUS  based 
quantitative  methodology  is  necessary  to  evaluate  both  the  final  system  perfor¬ 
mance  and  clinical  outcomes.  Research  is  currently  underway  to  remove  these 
limitations,  and  to  conduct  a  more  detailed  study  using  clinical  indicators. 
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Abstract 

Prostate  Brachytherapy  has  emerged  as  a  common  treatment  modality  for  early  stage  prostate  cancer, 
wherein  small  seeds  are  implanted  into  the  prostate  to  eradicate  the  cancer  by  emitting  radiation.  The 
main  limitation  in  contemporary  brachytherapy  is  faulty  seed  placement,  predominantly  due  to  the 
presence  of  intra-operative  edema  (tissue  expansion).  Intra-operative  measurement  of  edema  in  prostate 
brachytherapy  requires  localization  of  inserted  radioactive  seeds  relative  to  the  prostate.  Seeds  were 
reconstmcted  using  a  typical  non-isocentric  C-arm,  and  exported  to  the  commercial  brachytherapy 
delivery  system.  Technical  obstacles  for  3D  reconstmction  on  a  non-isocentric  C-arm  included  pose- 
dependent  C-arm  calibration;  distortion  correction;  pose  estimation  of  C-arm  images;  seed  reconstmction; 
and  C-arm  to  TRUS  registration. 

In  precision-machined  hard  phantoms  with  40-100  seeds  and  soft  tissue  phantoms  with  45-87  seeds,  we 
correctly  reconstmcted  the  seed  implant  shape  with  an  average  3D  precision  of  0.35  mm  and  0.24  mm, 
respectively.  The  reconstmcted  3D  seed  positions  were  then  registered  to  the  pre-plan  computed  from 
TRUS.  In  a  Phase- 1  clinical  trial  on  6  patients  with  49-82  planned  seeds,  we  achieved  intra-operative 
monitoring  of  seed  distribution  and  dosimetry.  We  optimized  the  100%  prescribed  iso-dose  contour  by 
inserting  an  average  of  4.17  additional  seeds,  making  intra-operative  detection  of  edema  (and  subsequent 
dose  optimization)  possible  on  a  typical  non-isocentric  C-arm,  at  negligible  additional  cost  to  the  existing 
clinical  installation.  Furthermore,  the  proposed  system  is  the  first  of  its  kind  that  can  intra-operatively 
detect  seed  migration  patterns.  For  many  seeds,  the  results  indicate  a  migration  greater  than  10  mm  from 
their  intended  dropped  locations. 

Keywords:  Prostate  Brachytherapy,  Edema,  Intra-operative  Dosimetry. 

Description  of  purpose:  Prostate  cancer  continues  to  be  a  significant  health  problem,  both 
nationally  and  worldwide.  Permanent  low  dose  rate  brachytherapy  has  emerged  as  a  definitive 
treatment  modality  of  early  stage  low  risk  prostate  cancer.  The  main  limitation  in  contemporary 
brachytherapy  is  faulty  seed  placement,  which  may  cause  insufficient  dose  to  the  cancer  cells 
and/or  inadvertent  radiation  of  the  rectum,  urethra,  and  bladder.  One  of  the  main  reasons  is 
edema,  a  natural  response  of  the  body,  wherein  the  prostate  tissue  expands  during  the  surgery  owing  to 
repeated  needle-insertions  and  high  radiation  dose  delivered  to  the  tissue.  It  is  one  of  the  leading  causes 
for  under-dosed  regions  in  the  prostate,  rendering  any  pre-operative  planning  imperfect.  LDR  prostate 
brachytherapy  is  almost  exclusively  performed  with  transrectal  ultrasound  (TRUS)  guidance, 
while  C-arm  fluoroscopy  is  ubiquitously  applied  for  gross  visual  assessment  of  the  implant.  C- 
arms  arms  are  ubiquitous  in  contemporary  LDR  prostate  brachytherapy,  yet  they  cannot  be  used 
for  real-time  intra-operative  dosimetric  analysis  due  to  a  plethora  of  technical  limitations.  We 
developed  a  system  that  overcame  these  limitations  toward  providing  quantitative  dosimetry 


analysis  and  implant  optimization.  In  this  paper  we  report  novel  method  and  system  to 
reconstruct  and  spatially  register  the  implanted  seeds  (that  are  visible  in  fluoroscopy)  and  soft 
tissue  anatomy  (that  is  visible  in  TRUS). 

Method:  We  employ  a  regular  clinical  brachytherapy  setup,  without  alteration,  including  a 
treatment  planning  workstation  and 
(B&K  Medical),  C-arm  (GE  OEC, 

Salt  Lake  City,  UT).  The  C-arm  is 
interfaced  with  a  standalone  laptop 
computer  through  an  NTSC  video 
line  and  frame  grabber. 

Workflow:  The  clinical  workflow  is 
identical  to  the  standard  procedure 
until  the  clinician  decides  to  run  a 
reconstruction  and  optimization, 
typically  at  each  half  of  the 
procedure.  A  set  of  C-arm  images 
are  collected  with  as  wide  separation 
as  possible  (15°)  and  transferred  to 
the  laptop.  After  processing  the 
fluoro  images,  the  seeds  are 
reconstructed  in  TRUS  space  and 
exported  to  the  Interplant®.  The 
physician  uses  standard  Interplant®  features  to  analyze  the  dose  and  modify  the  remainder  of  the 
implant  plan,  in  order  to  optimize  it.  The  procedure  concludes  when  the  exit  dosimetry  shows  no 
cold  spots. 


Numerous  technical  obstacles  had 
to  be  overcome  in  realizing  the 
workflow  above:  (a)  C-arm 
calibration;  (b)  image  distortion 
correction;  (c)  pose  estimation  of 
C-arm  images;  (d)  registration  of 
C-arm  to  TRUS;  (e)  seed 
reconstruction,  (f)  dosimetry 
analysis,  and  finally  (g)  implant 
optimization. 

Experiments  and  Results:  We 

present  results  from  6  patients 
enrolled  in  a  Phase-I  clinical  trial. 

Implant  reconstruction  and  dose 
analysis  was  performed  halfway  of 
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Figure  2:  Intra-operative  visualization  of  the  edema  is  possible  with 
the  system.  The  4 planned ’  (green)  versus  the  4 reconstructed ’  (red)  seed 
positions  as  seen  in  the  X-Y  plane  (template  view).  The  image 
acquisition  is  at  end  of  procedure.  A  trend  of  outward  radiation  of  the 
seeds  from  their  initial  seed  position  is  observed.  This  is  due  to  edema 
occurring  during  the  procedure. 
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Figure  1:  Overview  of  the  proposed  solution.  The  FTRAC 
fiducial  tracks  C-arms,  and  also  registers  TRUS  to  C-arm  images, 
making  quantitative  brachvtherapv  possible  on  anv  C-arm. 


the  implant  and  also  in  the  end  (an  additional  one  was  used  for  a  research  comparison  to  post-op 
CT).  In  all  the  patients,  the  final  dosimetry  detected  one  or  more  cold  spots  (Figure  3).  In  the 
fourth  patient,  9  additional  seeds  (initial 
67)  were  inserted  to  bring  up  VI 00  dose  to 
the  desired  level.  Apparently,  the 
physician  grew  quickly  to  trust  the  system 
in  catching  cold  spots,  and  instead  became 
very  cautions  not  to  create  hot  spots. 

Importantly,  all  patients  were  released 
from  the  operating  room  with  satisfactory 
V  I 00/V90  dose  and  favorable  DVH  for  all 
relevant  structures.  The  use  of  the  system 
added  about  45  minutes  to  the  OR  time, 
with  each  set  of  three  implant 
reconstructions  taking  about  15  minutes, 
including  image  capture,  processing,  3D 
reconstruction,  and  dosimetry  evaluation. 

We  anticipate  needing  only  a  single 
reconstruction  in  the  post-research  phase. 

Intra-operative  visualization  of  edema  was 
also  possible  with  the  system,  which  was 
found  to  be  significant  (Figure  2).  The 
seeds  showed  a  clear  tendency  of  outward 
migration  from  their  drop  positions  as  the 
procedure  progressed.  Furthermore,  post¬ 
operative  edema  can  also  be  computed 
from  a  comparison  of  intra-op  seed 
coordinates  to  their  post-implant  CT  coordinates. 

Conclusions:  The  system  was  exquisitely  accurate  in  technical  trials  and  has  shown  usefulness 
and  great  potential  in  a  limited  Phase- 1  trial.  It  can  quantitatively  measure  the  amount  of  edema 
and  the  associated  patterns  observed  in  a  patient. 

New  or  breakthrough  work  to  be  presented:  A  first  of  a  kind  system  that  can  use  a  non-isocentric  C- 
arm  to  intra-operatively  measure  the  pattern  and  extent  of  seed  migration  in  prostate  brachytherapy. 

Indicate  whether  the  work  is  being,  or  has  been,  submitted  for  publication  or  presentation 
elsewhere,  and,  if  so,  indicate  how  the  submissions  differ:  Some  preliminary  results  of  the  whole 
system  have  been  accepted  for  publication  at  MICCAI  2007.  This  submission  includes  complete  and 
detailed  results  with  extensive  analysis  from  all  the  patients  enrolled  in  the  Phase-I  clinical  trial. 


Figure  3:  The  system  is  able  to  detect  cold  spots.  The 
VI 00  contours  (blue)  as  assumed  by  the  planning  system 
(top)  and  as  computed  by  the  proposed  system  (bottom).  In 
this  patient,  3  cold  spots  were  discovered  in  exit  dosimetry. 
The  prostate  boundary  is  delineated  by  dark  red  contours 
and  the  green  squares  mark  current  seed  coordinates.  Note 
that  the  system  also  detected  4  seeds  that  have  drifted  out 
of  slice,  owing  largely  to  the  swelling  effects  of  edema. 
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Abstract.  We  developed  a  discrete  tomography  method  for  prostate  implant 
reconstructions  using  only  a  limited  number  of  X-ray  projection  images.  A  3D 
voxel  volume  is  reconstructed  by  back-projection  and  using  distance  maps 
generated  from  the  projection  images.  The  true  seed  locations  are  extracted 
from  the  voxel  volume  while  false  positive  seeds  are  eliminated  using  a  novel 
optimal  geometry  coverage  model.  The  attractive  feature  of  our  method  is  that 
it  does  not  require  exact  seed  segmentation  of  the  X-ray  images  and  it  yields 
near  100%  correct  reconstruction  from  only  six  images  with  an  average 
reconstruction  accuracy  of  0.86  mm  (std=0.46mm). 


1.  Introduction 

Brachytherapy  is  a  definitive  treatment  for  low  risk  prostate  cancer  that  represents  the 
vast  majority  of  new  cases  diagnosed  nowadays.  The  brachytherapy  procedure  entails 
permanently  implanting  small  radioactive  seeds  into  the  prostate.  The  main  limitation 
of  contemporary  brachytherapy  is  faulty  seed  placement  that  may  result  in  insufficient 
dose  to  the  cancer  and/or  inadvertent  radiation  to  the  rectum,  urethra,  and  bladder. 
Intra-operative  implant  optimization  promises  a  major  clinical  breakthrough,  but  for 
this  technique  to  succeed  the  implanted  seeds  must  be  reconstructed  and  registered 
with  the  anatomy  [1].  This  work  concentrates  on  the  first  problem,  reconstruction. 

C-arm  X-ray  fluoroscopy  is  the  gold  standard  in  observing  brachytherapy  seeds 
and  therefore  is  a  natural  candidate  for  implant  localization.  The  3D  coordinates  of 
the  seeds  can  be  calculated  from  multiple  X-ray  images  upon  solving  the 
correspondence  problem  [2-6].  These  methods  uniformly  require  that  seeds  are 
accurately  segmented  in  the  X-ray  images.  Significant  research  has  been  dedicated  to 
this  problem,  still  without  clinically  robust  and  practical  solution.  To  make  the 
problem  worse,  typically  7%,  but  often  as  much  as  43%  of  the  seeds  can  be  hidden  in 
the  X-ray  images  [7],  and  the  recovery  of  these  seeds  is  an  exigent  task  that  often 
leaves  seeds  undetected.  Su  et  al.  [7]  proposed  a  solution  to  the  hidden  seed  problem 
by  resolving  seed  clusters  and  extending  previously  published  approach  [3],  but  it  still 
required  perfectly  localizing  all  visible  seeds  all  projection  images. 

Classic  tomosynthesis  might  seem  a  suitable  reconstruction,  but  unfortunately,  it  is 
impractical  in  brachytherapy  because  (1)  the  swing  space  of  the  C-arm  is  limited  due 
to  collision  hazards  and  (2)  the  number  of  X-ray  images  is  strictly  limited  in  order  to 
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save  the  patient  and  the  OR  crew  from  excessive  toxic  radiation.  Tutar  et  al  [8]  has 
proposed  a  variant  of  tomosynthesis  denoted  as  selective  back  projection.  However 
this  method  demands  a  large  number  of  images  (>7)  and  wide  C-arm  angle  (>25°)  to 
succeed.  It  is  also  prone  to  introducing  false  positive  (FP)  seeds,  which  from  a 
dosimetric  point  of  view  are  more  troublesome  than  hidden  seeds,  because  they  act 
toward  underdosing  the  cancer.  Tutar  et  al.  use  a  heuristic  rule  to  recognize  FP  seeds 
by  their  sizes,  but  since  C-arm  pose  estimation  and  calibration  errors  affect  the  size  of 
objects,  this  may  result  in  faulty 
separation  of  the  true  and  false 
seeds. 

Our  approach  using  discrete 
tomography  is  different  in  that 
after  generating  a  3D  volume 
using  back-projection  we  detect 
and  remove  the  false  positive 
seeds  by  solving  an  optimal 
coverage  problem.  We  achieved 
high  reconstruction  rate  with 
fewer  images. 


2.  Method 


Fig.  1.  Two  false  positive  seeds  are  introduced  by  three 
legitimate  seeds  when  their  projectiles  intersect  the  same  voxels. 


In  the  example  in  Figure  1,  three  X-ray  images  are  used.  Three  seeds  project  in  each 
image  and  leave  a  “mark”  in  the  image,  typically  a  dark  blotch  in  fluoroscopy.  After 
reconstructing  a  3D  voxel  volume  with  common  back-projection,  five  candidate  seeds 
are  found  in  the  volume  and  all  appear  to  be  legitimate  seeds  in  each  image.  The 
question  is  how  to  separate  the  true  seeds  from  the  false  positives.  In  this  simplistic 
example,  it  is  very  easy  to  identify  the  legitimate  seeds  (solid  circles)  because  any 
other  choice  will  lead  to  an  inconsistency:  there  will  be  “seed  marks”  in  one  or  more 
X-ray  images  to  which  no  seed  in  the  volume  projects.  Therefore,  the  intuition  behind 
the  reasoning  is  that  each  “seed  mark”  in  each  image  must  be  “covered”  by  at  least 
one  seed  in  the  voxel  volume.  To  support  this  intuition,  we  develop  a  theoretical 
framework  based  on  optimal  geometry  coverage. 

Tomosynthesis  with  Distance  Map  and  Seed  Localization.  Tomosynthesis 
[10]  is  the  technique  to  reconstruct  a  3D  volume  from  multiple  2D  projection  images 
within  a  limited  angle.  The  most  commonly  used  approach  is  the  back  projection 
method,  which  was  first  introduced  for  CT  reconstruction.  In  this  method,  each  voxel 
in  a  3D  volume  is  projected  onto  all  the  images.  The  value  assigned  to  this  voxel  is 
calculated  as  the  average  of  the  intensity  values  of  its  projected  locations  in  the 
images. 

A  local  coordinate  system  must  be  defined  for  the  C-arm.  For  any  arbitrary  3D 
point  in  the  space,  after  C-arm  calibration,  its  projected  coordinates  on  the  2D  image 
plane  can  be  calculated  by  the  rules  of  perspective  projection. 

In  practice,  the  C-arm’ s  pose  is  estimated  with  some  error.  Since  the  size  of  the 
seeds  is  small  relative  to  the  focal  length,  the  reconstruction  is  quite  sensitive  to  pose 
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error.  To  make  the  reconstruction  more  robust  to  pose  error,  in  the  reconstruction  of 
the  3D  voxel  volume  we  use  distance  maps  rather  than  the  projection  image  itself. 

As  preprocessing,  in  each  image  we  first  extract  the  so  called  2D  seed  regions— 
areas  that  contain  seeds’  projections— using  adaptive  thresholding  and  morphological 
operators  and  call  the  resulting  images  “seed-only”  image.  Then  for  each  seed-only 
image,  a  distance  map  is  calculated  using  a  distance  transform:  the  value  at  each  pixel 
is  the  Euclidean  distance  to  its  nearest  2D  seed  region.  (Pixels  inside  a  2D  seed  region 
all  take  a  value  of  zero.)  In  reconstructing  the  3D  volume,  the  value  of  a  voxel  is  the 
average  of  the  distance  values  at  all  of  its  projected  locations. 

After  the  3D  voxel  volume  is  reconstructed,  candidate  3D  seed  regions  are 
extracted  by  thresholding.  The  threshold  value  is  based  on  estimated  pose  error.  For 
example,  in  case  of  small  pose  error  of  less  than  1  degree  rotation,  such  as  in  Jain  et 
al.  [6],  the  threshold  is  set  to  lA  pixel.  Upon  thresholding,  connected  3D  seed  regions 
are  considered  as  candidate  seeds.  The  candidate  seeds  are  then  labeled  using  the 
standard  3D  “connected  component  labeling”  method,  i.e.  any  two  neighboring  voxels 
in  the  3D  seed  regions  are  assigned  the  same  label  and  considered  as  parts  of  the  same 
seed.  After  that,  the  centroids  of  all  the  candidate  seeds  are  calculated  by  averaging 
the  3D  coordinates  of  all  voxels  wearing  the  same  label.  Each  candidate  seed  is  then 
represented  by  the  location  of  its  centroid,  in  addition  to  its  label. 

As  we  mentioned  earlier,  a  decisive  advantage  of  tomosynthesis  over  the  three-film 
technique  and  its  derivatives  [6]  is  that,  in  addition  to  not  requiring  precise  seed 
segmentation,  it  can  reconstruct  all  hidden  seeds.  But  the  disadvantage  is  that 
tomosynthesis  introduces  false  positive  seeds,  often  as  much  as  20%.  Next,  we 
propose  a  theoretical  framework  for  separating  the  true  seeds  from  the  false  ones. 

Theoretical  Framework  for  Separating  True  Seeds.  The  problem  is 
formulated  as  an  optimal  geometric  coverage  problem  [11].  The  general  optimal 
coverage  problem  arises,  for  example,  in  wireless  network  design.  From  a  given  set  of 
“server”  points,  we  must  select  the  minimum  subset  that  can  cover  a  given  set  of 
“client”  points.  The  goal  is  formulated  as  to  minimize  a  total  cost  function,  which  is 
the  sum  of  the  cost  functions  defined  on  all  the  selected  server  points. 

In  our  problem,  based  on  the  intuition  mentioned  earlier,  we  want  to  find  the  M 
true  seeds  from  the  set  of  N  candidates,  such  that  all  the  2D  seed  regions  are 
“covered”  in  all  projection  images,  i.e.  every  2D  seed  region  must  have  at  least  one 
true  seed  projected  in  it.  In  constructing  the  appropriate  cost  function,  we  can  use  the 
observation  that  a  false  positive  (FP)  seed,  owing  to  its  very  nature,  is  projected  close 
to  some  true  seed  in  every  image.  While  a  true  seed  may  be  projected  close  to  some 
other  true  seeds  in  some  images,  it  usually  is  not  the  case  for  every  image.  (While  it  is 
common  that  a  true  seed  can  be  hidden  in  some  projection  images,  we  have  never 
encountered  a  situation  where  a  true  seed  was  hidden  in  all  images.  We  also  note  that 
no  existing  method  can  recover  such  a  seed.) 

Hence  for  the  server  model,  we  define  the  cost  function  of  a  given  seed  (potential 
server  point)  as  the  sum  of  the  closest  distances  between  the  projections  of  this  seed 
and  the  projections  of  all  other  true  seeds,  in  all  images.  We  formulate  the  problem  as 
such  a  general  optimal  coverage  problem:  Given  the  N  candidate  seeds,  we  want  to 
find  M  seeds  ( x1 ,  x2 ,  . . .,  xM )  from  them  such  that  all  2D  seed  regions  in  the  seed- 
only  images  are  covered,  and  the  below  cost  function  is  minimized 
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M  n  A 

C  =  -1t'Ldi(xm)  (2) 

m=\i=\ 

where  di(xn)  =  min  ||  PJxm -Pjxn  ||,  V«  e  [1,2,--*,M], and/?  *  m  (3) 

m 

and  Pj  is  the  projection  operator  that  projects  a  3D  point  onto  the  jth  image. 

Unfortunately,  the  optimal  coverage  problem  is  NP-hard  [12]  and  its  computational 
complexity  is  0( C £ ),  where  C means  N  choose  M.  Thus  for  a  large  set  of  seeds,  it 
is  not  possible  to  find  its  global  optimal  solution.  We,  however,  managed  to  reduce 
the  size  of  the  problem  by  using  the  2D  seed-only  images  for  regularization.  To 
further  reduce  the  problem,  we  also  used  greedy  search  to  minimize  local  costs  rather 
than  the  global  cost  function  in  Eq.  (2). 


Seed  Clustering.  To  ensure  that  the  projections  of  selected  true  seeds  cover  all  the 
2D  seed  regions  in  all  seed-only  images,  the  candidate  seeds  are  clustered  based  on 
their  projections  in  each  image.  For  this  purpose,  we  label  all  2D  seed  regions  in  all 
seed-only  images,  in  the  similar  way  as  3D  labeling  described,  e.g.  the  separated  seed 
regions  are  assigned  different  labels  and  pixels  in  the  same  connected  region  wear  the 
same  label.  The  projections  of  seeds  can  then  be  clustered  based  on  these  labels.  An 
example  is  shown  in  Fig.  2.  Unlike  Su  et  al.  [7],  the  purpose  of  seed  clustering  is  not 
to  segment  the  2D  seed  projections  from  2D  seed  clusters  in  the  images,  or  to  identify 
the  number  of  true  seeds  in  each  cluster.  Instead,  we  use  the  seed  clustering  as  a  way 
to  relating  the  3D  seeds,  and  the  relationship  is  used  in  the  coverage  function  that  is  to 
be  minimized. 


Fig.  2.  (Left)  X-ray  image  on  phantom  data;  (Middle)  the  seed-only  image  resulted  after 
preprocessing;  (Right)  Example  of  seed  region  labeling  and  seed  grouping. 


On  a  projection  image,  the  projections  of  all  candidate  seeds  are  first  computed 
using  Eq.  (1).  We  denote  a  3D  candidate  seed  as  x„}  and  its  2D  projection  on  the  jth 
image  is  Pjxn.  For  each  Pjxn,  we  find  its  nearest  seed  region,  which  is  labeled 
as Lj(Pjxn) .  The  distance  fromE^  to  its  closest  seed  region  is  also  calculated  and 
denoted  as  dJn  ( d{  =  0  if  Pjxn  is  inside  a  seed  region).  The  seed  projections  are  then 
clustered  based  on LJ  (Pjxn) .  Let  there  be  Kj  seed  regions  on  the  jth  seed-only  image, 
and  they  are  labeled  as  1,  2,  KJ .  The  projections  of  the  N  candidate  seeds  on  the 
jth  image  are  then  clustered  into  Kj  sets  £l{  ,  such  that 

\/Pl,p2e  £li,L\Px)  =  Lj(Px\  j  =  1,  2,  ...,  Np,  k  =  l,2,...,K^  (4) 

Let  ||  ||  be  the  cardinal  of  setQ^  ,  we  say  the  seed  region  with  label  k  in  jth  image  is 
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covered  by  ||  ||  seeds.  \\ClJk  ||>  1 .  The  clustering  is  repeated  on  all  images. 

In  the  example  in  Fig.  2(c),  the  seed  regions  are  labeled  from  1  to  5  and  red 
asterisks  mark  the  2D  projections  of  the  candidate  seeds.  The  three  regions  labeled 
with  1,  4  and  5  contain  one  true  seed  each.  The  two  regions  labeled  with  2  and 
3 contain  two  (or  more)  true  seeds  each.  Since  the  seeds  in  region  2  (and  3)  are 
connected,  it  is  considered  as  a  single  seed  region.  The  candidate  seeds  include  both 
true  seeds  and  a  few  false  seeds.  In  this  example,  upon  clustering,  the  sets  with  label 
from  1  to  5  have  2,  3,  4,  1,  and  1  element(s)  individually 

Local  Coverage  and  Greedy  Search.  The  seed  clustering  according  to  each 
projection  can  help  reduce  the  size  of  our  optimization  problem.  If  a  seed  region  is 
covered  only  by  one  seed,  this  seed  must  be  a  true  one,  because  otherwise  this  region 
is  not  covered.  The  set  of  such  seeds  can  be  expressed  as: 

G  =  \Jj{xn  :  L' (P'xn)  =  k  and|| QJk  ||=1,  k  =  1,2,...,*'},  j  =  lZ...,Np  (5) 

where  Np  is  the  number  of  projection  images.  These  seeds  are  always  chosen  as  true 
seeds.  Hence  the  optimization  problem  is  reduced  to  choose  (M-\\G\ \)  true  seeds  from 
(7V-|  |G|  \)  candidate  seeds.  This  can  also  be  seen  in  the  example  in  Fig.  2(c).  The  two 
regions  labeled  with  4  and  5  contain  only  one  seed  each.  Therefore  these  two  seeds 
are  considered  as  true  seeds  and  are  always  chosen. 

Besides,  instead  of  finding  the  global  optimization  of  the  cost  function  in  (2),  we 
use  greedy  search  to  find  an  approximate  optimal  solution.  We  also  redefine  the  local 
coverage  cost  function  as 

NP  i  +  nJ  (x  1 

C  =  Tc(xn\  c{xn)  =  -YJ - — ,  forallx„£G  (6) 

n  7=1  1  +  dJn 

where 

DJ  ( xn )  =  min  ||  PJxn  -  PJxm  ||,  Mxm  eS-G,  and  V  ( Pjx„ )  =  V  ( Pjxm )  (7) 

m^n 

is  the  minimum  distance  between  the  projection  of  xn  and  the  projections  of  other 
seeds  in  the  cluster  in  the  jth  image  that  includes  xn,  and  dJn  is  the  distance  from  Pjxn 

to  the  nearest  seed  region.  dJn  is  added  in  the  cost  function  to  include  the  effect  of 
imperfect  X-ray  pose  estimation. 

The  minimization  problem  in  (6)  is  solved  using  greedy  search  iteratively.  During 
each  cycle  of  iteration,  the  seed  that  has  the  largest  cost  value  c(xn )  is  considered  as  a 
false  seed  and  is  removed  from  the  candidate  seed  set.  G  is  updated  at  the  beginning 
of  each  cycle  of  iteration,  for  after  the  removal  of  one  seed,  there  may  be  additional 
seeds  that  cover  some  region  alone  (i.e.  no  other  seed  covers  this  region.)  These  seeds 
need  to  be  extracted  and  added  to  G.  The  algorithm  is  summarized  as  below. 
Algorithm  1 :  find  M  good  seeds  from  N  candidate  seeds  using  greedy  search 

1 .  Initialize  S  be  the  set  of  candidate  seeds. 

2.  For  /  =  1:N-M 

3.  Calculate  G  using  (5). 

4.  Calculate  local  cost  function  c(x)  of  all  seeds  x  in  S  -  G  using  (6). 

5.  Find  xk  e  S  -  G  ,  such  that  xk  =  argx  min  c(x) . 

6.  Remove  xk  from  S:  S  =  S  -  {x^} . 
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3.  Experiments  and  Results 


Fig.  3.  (Left)  Part  of  a  synthetic  C-arm  image;  (Right)  One  slice  of  the  reconstructed  3D 
volume.  The  red  circles  mark  the  true  seeds,  and  the  green  squares  represent  the  FP  seeds. 

Simulations.  Synthetic  C-arm  images  were  used  to  verify  our  method.  The  images 
simulated  a  50  cc  prostate  with  a  seed  density  of  2.0  seed/cc.  The  C-arm’ s  focal 
length  was  1000  mm,  and  the  pixel  size  was  0.25  mm.  Six  images  were  generated  on 
a  20°  cone  around  the  AP  axis  with  evenly  distributed  angles.  The  seeds  were 
represented  by  cylinders  with  a  radius  of  0.4  mm  and  a  length  of  1 .45  mm.  A  typical 
synthetic  image  is  shown  in  Fig.  3(a).  No  pose  estimation  error  was  assumed. 


Table  1.  Seed  reconstruction  results  using  different  number  of  synthetic  images 


Number  of 
images 
used 

Number  of 
true  seeds 
implanted 

Number  of 
candidate  seeds 
before  FP  removal 

Correctly 

reconstructed  seeds 
after  FP  removal  (%) 

Mean 

reconstruction 
error  (mm) 

Reconstructs 
n  error  STD 
(mm) 

3 

96 

105.6 

99.4 

0.19 

0.19 

4 

96 

97.7 

100 

0.12 

0.09 

We  reconstructed  the  seeds  using  3,  4  and  6  images.  For  3  images,  20  experiments 
were  performed  by  using  all  the  20  combinations  of  selecting  3  images  from  the  6 
available  images.  For  4  images,  four  experiments  were  performed  using  different 
image  combinations.  The  results  are  shown  in  Table  1.  In  each  experiment,  the 
number  of  implanted  seeds  was  assumed  known.  After  FP  removal,  a  set  of  candidate 
seeds  that  equal  to  the  number  of  implanted  seeds  were  chose.  These  chosen  seeds 
were  then  compared  with  the  known  3D  locations  of  implanted  seeds.  When  using  3 
images,  there  were  about  10%  FP  seeds  reconstructed  initially  and  then  almost  all 
were  successfully  removed,  indicated  by  the  near  100%  correct  final  reconstruction 
rate.  This  means  that  when  there  is  no  pose  estimation  error,  the  seeds  can  be 
reconstructed  accurately  even  from  as  few  as  3  images  using  our  method. 

Phantom  Studies.  Experiments  were  performed  on  a  seed  phantom  constructed 
from  acetol.  The  FTRAC  fiducial  was  used  to  track  the  C-arm  (with  accuracy  of  0.56 
mm  translation  and  0.33°  rotation).  It  was  attached  to  the  seed  phantom,  as  shown  in 
Fig.  4.  The  seed  phantom  is  comprised  of  twelve  slabs  with  5mm  thickness  each. 
Each  slab  has  more  than  100  holes  with  5  mm  spacing,  into  which  the  seed  can  be 
inserted.  By  precise  manufacturing,  the  seed  phantom  was  attached  in  a  known 
position  to  a  radiographic  C-arm  tracking  fiducial,  replicated  after  Jain  et  al.  [13].  In 
this  way,  the  exact  location  of  each  implanted  seed  was  known  relative  to  the  tracking 
fiducial,  serving  as  ground  truth.  Because  some  rotation  error  got  introduced  during 
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assembly,  the  ground  absolute  seed  positions  had  about  0.5  mm  error.  (Note  that  their 
relative  positions  were  precisely  known.)  The  seed  density  was  about  1.56  seed/cc. 
Five  data  sets  were  collected  with  the  number  of  seeds  varying  from  40  to.  For  each 
data  set,  six  images  within  a  20°  cone  around  the  AP  axis  were  taken  using  a  Phillips 
integris  V3000  fluoroscope  and  dewarped  using  the  pin-cushion  test.  We  used  four, 
five  and  six  images  to  reconstruct  the  seeds,  by  using  all  the  combinations  of  the  six 
available  images,  i.e.  15  experiments  using  four  images,  6  using  five  images,  and  1 
using  six  images  for  each  data  set.  The  reconstructed  seeds  were  compared  with  the 
computed  ground  truth,  and  the  results  were  shown  in  Table  2. 


Fig.  4  (Left)  An  image  of  the  seed  phantom  attached  to  the  FTRAC  fiducial.  (Right)  A  typical 
X-ray  image  of  the  combination  (with  100  seeds  inserted) 


Table  2.  Seed  reconstruction  results  on  phantom  data  using  different  number  of  images 


Number  of 
images 
used 

Number  of 
seeds 
implanted 

Number  of 
candidate  seeds 
before  FP  removal 

Correctly 

reconstructed  seeds 
after  FP  removal  (%) 

Mean 

reconstruction 
error  (mm) 

Reconstruction 
error  STD 
(mm) 

6 

40 

46 

100 

1.04 

0.56 

55 

58 

100 

0.67 

0.39 

70 

82 

98.6 

0.72 

0.42 

85 

94 

100 

0.97 

0.44 

100 

112 

98 

0.94 

0.52 

5 

40 

46.3 

99.5 

1.02 

0.53 

55 

68.3 

98.8 

0.85 

0.46 

70 

92.3 

97.8 

0.85 

0.46 

85 

105.3 

97.0 

1.00 

0.57 

100 

121.3 

92.5 

1.19 

0.81 

4 

40 

53.7 

96.8 

1.24 

0.76 

55 

82.1 

94.7 

0.99 

0.69 

70 

112 

94.3 

0.91 

0.77 

85 

135.8 

90.1 

1.13 

0.62 

100 

159.7 

86.3 

1.44 

0.91 

It  was  shown  in  [12]  that  in  prostate  brachytherapy,  95%  of  the  implanted  seeds 
must  be  recovered  in  order  to  obtain  clinically  accurate  estimation  of  the  dose.  Our 
results  showed  that  (1)  from  six  images  more  than  98%  seeds  can  be  successfully 
reconstructed,  and  (2)  six  images  were  sufficient  for  deriving  clinically  accurate  dose 
estimation.  Furthermore,  the  number  of  required  images  depends  on  the  number  of 
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seeds  implanted.  For  a  smaller  number  of  seeds  (40,  55,  70),  only  four  to  five  images 
were  required.  Additionally  the  performance  of  our  method  depends  on  the  ratio 
between  the  number  of  FP  seeds  and  the  implanted  seeds,  which  in  turn  depends  on 
both  the  number  of  images  used  and  the  number  of  implanted  seeds.  Generally  the 
lower  is  this  ratio,  the  higher  is  the  seed  reconstruction  ratio  we  may  achieve. 

In  summary,  we  presented  a  novel  method  for  prostate  brachytherapy  seed 
reconstruction  using  C-arm  images.  We  generate  distance  maps  from  the  2D 
projection  images,  then  a  3D  volume  is  then  reconstructed  using  tomosynthesis  using 
the  distance  maps,  and  finally  true  seeds  are  extracted  from  the  voxel  volume.  The 
attractive  feature  of  our  method  is  that  it  does  not  require  exact  seed  segmentation  of 
the  X-ray  images.  As  a  tradeoff,  our  method  requires  slightly  higher  number  of 
images  than  the  methods  that  requires  elaborate  explicit  segmentation.  Our  method 
yields  near  100%  correct  reconstruction  from  only  six  images  with  an  average 
reconstruction  accuracy  of  0.86  mm  (std=0.46mm).  The  method  was  robust  to  pose 
error  present  in  radiographic  C-arm  tracking. 
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Abstract.  X-ray  C-arm  fluoroscopy  is  a  natural  choice  for  intra-operative 
seed  localization  in  prostate  brachytherapy.  Resolving  the  correspon¬ 
dence  of  seeds  in  the  projection  images  can  be  modeled  as  an  assignment 
problem  that  is  NP-hard.  Our  approach  rests  on  the  practical  observation 
that  the  optimal  solution  has  almost  zero  cost  if  the  pose  of  the  C-arm 
is  known  accurately.  This  allowed  us  to  to  derive  an  equivalent  problem 
of  reduced  dimensionality  that,  with  linear  programming,  can  be  solved 
efficiently  in  polynomial  time.  Additionally,  our  method  demonstrates 
significantly  increased  robustness  to  C-arm  pose  errors  when  compared 
to  the  prior  art.  Because  under  actual  clinical  circumstances  it  is  exceed¬ 
ingly  difficult  to  track  the  C-arm,  easing  on  this  constraint  has  additional 
practical  utility. 


1  Introduction 

Intraoperative  dosimetric  quality  assurance  in  prostate  brachytherapy  critically 
depends  on  discerning  the  three-dimensional  (3D)  locations  of  implanted  seeds. 
The  ability  to  reconstruct  the  implanted  seeds  intraoperatively  will  allow  us  to 
make  immediate  provisions  for  dosimetric  deviations  from  the  optimal  implant 
plan.  A  method  for  seed  reconstruction  from  pre-segmented  C-arm  fluoroscopy 
images  has  been  proposed,  among  other  works,  by  Jain  et  al.  in  [1],  where  the  3D 
coordinates  of  the  implanted  seeds  are  calculated  upon  resolving  the  matching 
of  seeds  in  multiple  X-ray  images. 

At  least  three  images  are  necessary  to  eliminate  ambiguities.  The  resulting 
optimization  problem  is  NP-hard.  Heuristic  approaches,  such  as  of  Jain  et  al.  [1], 
have  been  proposed  to  approximately  solve  the  optimization  problem.  However, 
the  use  of  a  heuristics  leads  to  algorithmic  error ,  in  addition  to  physical  errors 
like  the  inaccuracy  in  knowing  the  relative  poses  of  the  C-arm  shots  (pose  error). 
To  tackle  this  issue,  we  propose  to  consider  all  the  images  simultaneously,  instead 
of  suboptimal  subsets  of  two  images  such  as  proposed  in  [1]. 

The  optimization  problem  has  a  salient  feature:  since  the  images  represent 
a  real  situation  (i.e.  an  existing  object,  the  set  of  seeds,  is  being  imaged),  the 
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optimal  solution  has  a  near-zero  cost  when  the  pose  error  is  low,  and  is  actu¬ 
ally  zero  without  pose  error.  We  propose  to  utilize  this  feature  of  the  problem 
to  reduce  its  number  of  variables,  thereby  allowing  to  obtain  the  optimal  solu¬ 
tion  at  a  reasonable  computational  cost.  This  exact  dimensionality  reduction  is 
only  possible  when  the  pose  error  is  sufficiently  low.  We  claim  that  this  is  not 
restrictive  in  our  framework  since  a  high  pose  error  leads  to  high  error  in  the 
estimation  of  the  3D  coordinates  of  the  implanted  seeds,  which  is  not  acceptable. 
Actually,  the  idea  of  dimensionality  reduction  is  not  new.  For  instance,  in  [1,  p. 
3480]  the  original  tripartite  matching  is  projected  into  inspired  bipartite  match¬ 
ings,  while  introducing  inaccuracy.  In  contrast,  the  proposed  method  performs 
dimensionality  reduction  while  ensuring  equivalency  to  the  original  problem. 

The  MARSHAL  method  of  Jain  et  al.  has  demonstrated  solid  performance  [1] 
and  we  chose  this  work  as  the  benchmark  and  basis  of  comparison  for  our  work. 
A  comparison  between  MARSHAL  and  our  method  was  conducted  to  evaluate 
the  sensitivity  to  pose  errors  on  simulated  and  phantom  data.  The  proposed 
method  shows  significant  increase  in  robustness  to  pose  errors. 

2  Method 

Consider  a  collection  of  X-ray  images  of  a  constellation  of  implanted  seeds.  We 
assume  that  the  2D  seed  locations  can  be  identified  on  each  of  the  X-ray  images, 
and  we  consider  the  problem  of  identifying  corresponding  seeds  in  all  the  images. 
Given  these  matched  seed  locations,  a  reconstruction  of  the  seed  locations  in 
3D  can  be  achieved  provided  there  are  no  ambiguities.  It  is  more  likely  that 
such  ambiguities  are  avoided  when  there  are  more  X-ray  images,  but  this  in 
turn  increases  the  complexity  of  the  problem.  Here  we  do  not  consider  CT  or 
limited  angle  tomosynthesis,  as  our  work  focuses  on  reconstruction  from  a  very 
limited  number  of  images.  We  propose  a  solution  for  three  images,  which  is  often 
sufficient  in  practice,  and  which  is  extendable  to  more  images. 

2.1  3D  Reconstruction  as  a  matching  problem 

The  3D  locations  of  the  implanted  seeds,  modeled  as  points,  can  be  reconstructed 
through  3D  triangulation  from  the  X-ray  images  upon  resolving  the  correspon¬ 
dence  of  seeds,  which  is  the  focus  of  this  paper.  Let  n  be  the  number  of  points 
in  the  clinical  work  volume.  Let  be  the  position  of  Ith  point  in  mth  image. 
When  three  images  are  used,  the  matching  problem  can  be  formulated  as  an 
axial  3D  assignment  problem  (3DAP)  [1]: 

n  n  n 
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Cijk  is  the  the  cost  of  matching  point  sn  to  points  Sj2  and  Sk 3.  %ijk  is  a  binary 
variable  deciding  the  correctness  of  the  match  (i,  j,  k).  (2)  force  every  segmented 
point  to  be  chosen  once.  Thus,  x  represents  any  feasible  global  match,  with 
the  cost  of  that  correspondence  given  by  •  One  good  choice  for 

a  cost-metric  c  is  projection  error  (PE)  [1].  For  any  given  set  of  poses  and 
correspondence,  the  intersection  of  the  three  lines  that  join  each  projection  to 
its  respective  X-ray  source  can  be  computed  using  a  closed  form  solution  that 
minimizes  the  L 2  norm  of  the  error.  PE  can  be  computed  by  projecting  this 
3D  point  in  each  image  and  then  measuring  the  distance  between  the  projected 
location  and  the  observed  location  of  the  point. 

A  feasible  solution  x  of  the  above  problem  is  a  3D  permutation  array.  This 
problem  has  (n!)2  feasible  solutions.  Branch  and  bound  is  a  classical  algorithm 
for  optimally  solving  the  3DAP.  This  can  be  generally  achieve  only  for  n  small 
because  of  the  combinatorial  explosion.  Thus,  it  has  been  proposed  heuristics 
that  approximately  solve  the  3DAP,  such  as  MARSHAL  in  [1].  MARSHAL  sub- 
optimally  projects  the  original  3DAP  into  three  distinct  2DAP  that  can  be  solved 
in  polynomial  time  by  using  the  Hungarian  algorithm. 

We  point  out  that  the  3DAP  has  a  salient  feature  that  we  can  exploit.  Since 
the  images  represent  a  real  situation,  the  optimal  solution  has  a  near-zero  cost 
when  the  pose  error  is  low  and  the  optimal  cost  is  actually  zero  without  no  pose 
error.  In  the  next  section,  we  use  this  feature  to  reduce  the  number  of  variables 
in  the  problem,  thus  permitting  us  to  get  the  optimal  solution  at  a  reasonable 
computational  cost.  This  new  method  tackles  the  complete  optimization  problem 
without  using  suboptimal  projections,  such  as  in  MARSHAL. 


2.2  Dimensionality  reduction 

Let  N  =  n3.  We  rewrite  the  variables  Xijk  and  the  costs  in  a  vectorial  form 
such  that  x,  c  G  In  the  sequel  we  also  make  use  of  the  notation  U£  to  denote 
Uijk •  The  3DAP  (l)-(2)  reads  as  the  following  integer  linear  program 

P:  mine1  a?,  (3) 

xec  v  ' 

with  the  constraint  set  C  =  {x  :  M.x  =  [1, . . . ,  l]t ,  xg  E  {0, 1}},  where  Mx  = 
[1, . . .  ,1] 1  is  a  matrix  form  of  (2). 


Principle:  Since  the  coefficients  of  x  are  either  0  or  1  and  there  must  be  n 
l’s,  an  optimal  solution  of  P  can  be  thought  of  as  the  selection  of  n  cost  coeffi¬ 
cients  such  that  the  resulting  cost  is  minimum  while  constraints  C  are  satisfied. 
Given  a  feasible  solution,  Lemma  1  (below)  states  that  all  cost  coefficients  that 
are  greater  than  the  cost  associated  with  this  solution  cannot  be  selected  in  the 
optimal  solution.  Since  those  coefficients  can  never  be  selected,  the  dimension 
of  the  problem  can  be  reduced  by  removing  those  coefficients  from  further  con¬ 
sideration.  This  yields  an  equivalent  problem  of  reduced  dimensionality.  If  the 
reduction  in  dimensionality  is  sufficiently  large,  then  the  new  problem  can  be 


solved  exactly  in  reasonable  time  even  though  the  original  problem  is  far  too 
costly  to  solve. 


Lemma  1.  Let  us  assume  that  the  cost  coefficients  eg  are  positive.  Let  Xo  G  C  be 
a  feasible  solution.  The  integer  linear  problem  P  defined  by  (3)  is  equivalent  to 
the  following  integer  linear  problem  (ie.,  they  share  the  same  optimal  solutions) 

P'  :  mmfc'Yx,  where 
xecK  ' 

,  jce,  ifce<mP{x  0) 

ce  =  <  ^  ,  x  (4) 

l^oo,  ifce>mP(x  0) 

and  where  mp(x )  =  clx  is  the  cost  of  problem  P  at  the  feasible  solution  x. 
Proof.  Let  x*  G  C  be  an  optimal  solution  to  problem  P.  We  have  for  all  Xq  G  C 


mp(x*)  <  mp(x o).  (5) 

Let  us  consider  eg  >  mp{x o).  From  (5)  we  have  eg  >  mp{x*).  Since  eg  are 
positive  and  x%  are  binary,  we  have  necessarily  x\  =  0.  □ 


The  dimensionality  reduction  can  be  illustrated  as  follows: 


The  original  problem  P  is  equivalent  to  the  reduced  problem 


P  :  min  (fx, 
xec 


where  x,c  G  ( K  <  N )  and  with  the  constraint  set  C  =  {x  :  Mcc  = 
[1, . . . ,  l]1 ,  xg  G  {0, 1}}  with  M  =  MR  and  where  R  is  the  dimensionality  reduc¬ 
tion  matrix  of  size  N  x  K  such  that  [xn  0  x&  0  . . .  xik\  =  R  \x\  x^  •  •  •  xk\  • 
Once  the  reduced  problem  P  is  solved,  the  optimal  solution  to  the  original  prob¬ 
lem  P  is  simply  given  by  cc*  =  Rcc*. 

Note  that  a  dimensionality  reduction  ( K  <  N )  is  not  guaranteed  for  the 
general  3DAP,  even  in  the  most  favorable  case  mp(x o)  =  mp{x*).  However,  a 
dimensionality  reduction  occurs  in  the  case  of  our  problem  since  the  range  of  the 
cost  coefficients  is  wide  and  the  optimal  solution  has  a  near-zero  cost  when  the 
pose  error  is  low.  The  practical  interest  clearly  depends  on  the  dimensionality 
reduction  ratio  ( K/N ).  We  show  next  that  this  ratio  actually  can  be  improved. 


Improving  dimensionality  reduction:  Lemma  1  uses  only  the  integer  con¬ 
straint  in  C  to  reach  dimensionality  reduction.  But  using  the  fact  that  the  feasible 
set  C  comprises  permutation  matrices,  is  is  possible  to  reduce  the  dimensionality 
of  this  problem  even  further.  To  demonstrate  this,  we  start  with  the  following 
Lemma. 


Lemma  2.  The  integer  linear  problem  P  defined  by  (3)  is  equivalent  to  the 
following  integer  linear  problem 

P"  :  min(c//)tx 
xec 

where  the  minimum  cost  in  each  row  is  subtracted  from  the  entire  row 

c'lik=cijk-mmcijk.  (6) 

j,k 

Proof  For  lack  of  space,  the  proof  is  not  detailed. 

To  show  that  this  Lemma  permits  further  dimensionality  reduction,  let  us  apply 
Lemma  1  on  the  new  problem  P" .  From  (4),  the  cost  coefficient  c"  is  equivalent 
to  oo  if  c'f  >  mp"{x).  According  to  (6)  the  former  condition  is  equivalent  to 

n 

Cijk  >  mp(x )  -  (W  min  cijk  -  min cijk) 
z '  j,k  j,k 

i=  1 

The  latter  condition  is  less  restrictive  than  q  >  mp(x )  since  the  cost  coefficients 
are  positive.  As  a  result,  the  dimensionality  reduction  is  higher  within  the  new 
problem  P"  than  within  the  original  problem  P.  It  is  then  preferable  to  consider 
problem  P"  instead  of  problem  P  since  Lemma  2  ensures  that  they  are  equiva¬ 
lent.  It  is  actually  possible  to  reduce  dimensionality  even  further.  The  operation 
(6)  can  also  be  performed  successively  for  the  columns  and  depths  to  decrease 
the  cost  coefficients,  while  still  ensuring  equivalency  to  the  original  problem. 


2.3  Optimization  strategy 

Integer  Programming  /  Linear  Programming:  The  Integer  Program  (1)- 
(2)  can  be  directly  solved  with  standard  techniques  such  as  branch  and  bound. 
IP  problems  are  NP  hard,  however,  and  may  take  an  exponential  amount  of 
computational  time.  It  has  been  shown  that  the  linear  program  corresponding  to 
the  2D  assignment  problem  (2DAP)  has  an  integer  solution  even  without  integer 
constraints  [2].  As  well,  this  linear  program  can  be  solved  efficiently  in  polynomial 
time  using  interior  point  methods  for  instance  [3].  To  our  knowledge,  however, 
there  is  no  analogous  result  for  the  3DAP.  Nevertheless,  we  have  gone  ahead  and 
implemented  the  linear  program  for  the  3DAP  problem  followed  by  a  test  to  see 
if  its  solution  is  binary  (up  to  numerical  errors).  In  all  of  our  experiments,  we 
have  never  obtained  a  nonbinary  solution  to  this  problem,  which  points  to  the 
potential  validity  of  the  2D  theoretical  result  in  3D  as  well.  (We  are  currently 
investigating  this  theoretical  issue.) 

Dimensionality  reduction  thresholding:  From  Lemma  1,  the  degree  of  di¬ 
mensionality  reduction  depends  on  the  cost  mp(x o)  of  a  feasible  solution  xq.  It 
is  possible  to  use  a  suboptimal  algorithm,  such  as  MARSHAL,  to  find  a  feasible 
solution  xq.  Unfortunately,  when  there  is  high  pose  error,  MARSHAL  often  pro¬ 
vides  such  a  suboptimal  solution  that  very  little  dimensionality  reduction  can  be 


be  achieved.  Therefore,  we  propose  instead  to  choose  a  threshold  parameter  77, 
which  is  essentially  a  “guess”  as  to  what  the  cost  of  a  feasible  solution  might  be. 
This  permits  us  to  reduce  the  dimensionality  of  the  problem  and  run  the  linear 
problem  on  the  resultant  problem.  If  the  solution  of  this  problem  is  integer  and 
it  has  a  cost  lower  than  77,  then  it  must  be  optimal.  If  the  resultant  cost  is  larger 
than  77  then  the  solution  might  be  optimal,  but  we  cannot  guarantee  it.  It  is  then 
our  option  whether  to  accept  a  (potentially)  suboptimal  solution  or  to  increase 
77  and  rerun  the  linear  problem  until  we  have  a  guaranteed  optimal  solution. 

In  our  experiments,  we  determine  77  in  the  following  way.  Rank  order  all 
costs  from  lowest  to  highest  and  pick  an  integer  K.  Let  77  be  the  value  of  the 
Kth  smallest  cost  coefficient.  The  influence  of  K  on  the  rate  of  feasibility  and 
optimality  of  the  proposed  method  is  experimentally  studied  in  the  next  section. 

3  Simulation  and  Phantom  Experiments 

We  present  a  comparison  between  the  MARSHAL  algorithm  and  the  proposed 
method  using  simulation  and  phantom  data.  We  got  a  copy  of  the  MARSHAL 
code  from  the  authors  for  comparison  [1].  Both  algorithms  were  implemented  in 
Matlab  7.1  on  a  Linux  PC  (Intel  EM64T  3.6  GHz,  24GB  RAM). 


3.1  Evaluation  of  pose  errors  sensitivity 

Two  separate  comparisons  were  performed.  One  compared  the  two  algorithms 
to  translational  errors  and  the  other  to  rotational  errors.  These  are  both  com¬ 
mon  errors  in  C-arm  position  calibration  in  the  operating  room.  Random  error 
was  modeled  using  a  uniform  probability  density  function.  When  we  report  re¬ 
sults  for  an  h  error  this  means  that  each  of  the  three  components  of  error  (in 
either  translation  or  rotation)  were  generated  as  independent  random  variables 
uniformly  distributed  on  [—h,h\.  Following  [1],  a  statistical  bias  in  translation 
was  incorporated  in  the  generation  of  the  datasets  to  account  for  the  expected 
differences  in  directional  errors  in  fluoroscope  tracking  using  a  fiducial.  In  par¬ 
ticular,  we  assumed  that  the  in-plane  error  is  a  factor  of  five  times  smaller  than 
the  through-plane  error  h.  No  analogous  bias  was  used  in  the  rotation  errors. 

Realistic  simulations  of  prostate  brachytherapy  seeds  implants  were  gener¬ 
ated  with  seeds  density  of  2  and  2.5  seeds/cc  and  prostate  size  of  35  and  45  cc. 
The  number  of  seeds  in  the  implants  were  n  =  {72,  84,  96, 112}.  A  cone  angle  of 
10°  was  used  for  the  acquisition  of  the  three  simulated  X-ray  images.  Averaged 
results  from  a  total  of  2,  000  datasets  are  shown  in  Fig.  1. 

The  proposed  method,  with  lOOn  cost  coefficients  remaining  after  dimen¬ 
sionality  reduction,  performs  significantly  better  than  MARSHAL,  as  shown  in 
Fig.  l(a)-(d).  It  still  matches  correctly  89%  of  the  seeds  when  rotation  error 
reaches  4°,  while  MARSHAL  drops  to  59%.  The  proposed  method  still  matches 
correctly  99%  of  the  seeds  when  the  translation  error  reaches  10  mm,  while 
MARSHAL  drops  to  72%.  The  computational  time  of  the  proposed  method  is 
higher  than  that  of  MARSHAL  as  shown  in  Fig.  l(e)-(f).  We  point  out  that 
most  of  the  computational  time  of  the  proposed  method  (solid  line)  is  actually 
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Table  1.  Feasibility  and  optimality  of  the  proposed  method  for  pose  error. 


Feasibility  rate 

j  Guaranteed  optimality  rate 

Number  of  cost  coef. 

20n 

50n 

lOOn 

500n 

lOOOn 

20n 

50n 

lOOn 

500n 

lOOOn 

Error  1  ° 

0.9 

0.96 

1 

1 

0.98 

0.88 

0.91 

0.94 

0.98 

1 

Error  2  ° 

0.9 

0.98 

1 

0.98 

0.94 

0.35 

0.38 

0.46 

0.6 

0.71 

Error  3  ° 

0.58 

0.94 

0.98 

0.94 

0.92 

0.1 

0.13 

0.13 

0.22 

0.27 

Error  4  ° 

0.42 

0.83 

0.92 

0.92 

0.85 

0.15 

0.1 

0.07 

0.07 

0.1 

Error  2  mm 

0.96 

1 

1 

1 

0.94 

1 

1 

1 

1 

1 

Error  4  mm 

0.94 

0.98 

1 

1 

0.92 

0.56 

0.6 

0.  65 

0.77 

0.89 

Error  6  mm 

0.96 

0.96 

1 

0.96 

0.88 

0.35 

0.41 

0.46 

0.63 

0.64 

Error  8  mm 

0.77 

0.9 

0.98 

0.94 

0.88 

0.14 

0.19 

0.19 

0.31 

0.38 

Error  10  mm 

0.58 

0.92 

0.96 

0.94 

0.77 

0 

0.02 

0.04 

0.07 

0.08 

used  in  the  computation  of  the  cost  coefficients.  However,  computing  all  the  cost 
coefficients  is  not  required  since  most  of  them  will  eventually  be  thrown  out  by 
dimensionality  reduction.  After  further  code  optimization,  we  expect  the  com¬ 
putation  time  to  reduce,  near  to  the  time  required  solely  for  cost  minimization 
(dotted  line). 

The  feasibility  and  optimality  rates  of  the  proposed  method  as  a  function 
of  the  number  of  remaining  cost  coefficients  after  dimensionality  reduction  are 
shown  in  Tab.  1.  It  is  expected  that  the  feasibility  rate  should  increase  as  a 
function  of  the  number  of  cost  coefficients.  This  is  true  for  smaller  numbers  of 
cost  coefficients  but,  surprisingly,  the  feasibility  rate  decreases  when  the  number 
of  cost  coefficients  reaches  500n  (~  50,  000).  This  is  due  to  the  failure  of  linprog 
in  Matlab  using  default  parameters  because  “one  or  more  of  the  residuals,  duality 
gap,  or  total  relative  error  has  stalled”.  These  cases  were  not  displayed  in  Fig.  1, 
and  we  are  currently  investigating  how  to  cope  with  them. 

The  guaranteed  optimality  rate  increases  as  a  function  of  the  number  of  cost 
coefficients.  For  low  errors  (1°  and  2  mm),  all  solutions  are  guaranteed  optimal 
given  an  acceptable  dimensionality  ratio.  We  point  out  that  the  guaranteed  opti¬ 
mality  is  only  a  sufficient  condition.  As  shown  in  Fig.  1(b),  all  solutions  from  0  to 
8  mm  translation  error  of  the  proposed  method  correspond  to  perfect  matching 
(100%),  even  if  they  are  not  all  guaranteed  optimal. 


3.2  Phantom  Experiments 

A  radiographic  fiducial  was  used  to  track  the  C-arm  (0.56  mm  translation;  0.33° 
rotation  accuracy)  and  was  accurately  attached  to  a  random  point  cloud  phan¬ 
tom.  Phantoms  of  {40,  55,  70,  85, 100}  points  with  1.56  points /cc  were  used.  Six 
X-ray  images  within  a  20°  cone  around  the  AP-axis  were  randomly  taken  using 
an  Philips  Integris  V3000  fluoroscope  and  dewarped.  Thus  both  the  seed  loca¬ 
tions  and  X-ray  pose  were  not  biased/optimized  in  any  way,  closely  representing 
an  uncontrolled  surgical  scenario.  Each  image  was  hand  segmented  to  establish 
the  true  segmentation  and  correspondence. 

Both  MARSHAL  and  the  proposed  method  perform  very  well  on  phantoms, 
achieving  almost  perfect  matching,  as  shown  in  Tab.  2.  Note  that  the  accuracy 
of  the  radiographic  fiducial  ensures  here  a  low  pose  error.  The  proposed  method 
is  significantly  slower  compared  to  MARSHAL.  We  point  out  that  the  proposed 
method  uses  here  more  than  90%  of  the  computational  time  for  the  n3  cost 
coefficients.  We  expect  to  reduce  significantly  this  time  as  explained  in  Sect.  3.1. 


Table  2.  Performance  of  MARSHAL  and  the  proposed  method  on  phantoms. 


MARSHAL  | 

|  Proposed  method 

Number  of  seeds 

40 

55 

70 

85 

100 

40 

55 

70 

85 

100 

Mean  Match.  (%) 

97.6 

100 

98 

97.7 

98.2 

98 

99.4 

97.1 

100 

98 

STD  Match.  (%) 

3.6 

0 

2.3 

3.2 

2.3 

2.6 

1.4 

0 

0 

0 

Time  (s) 

0.3 

0.6 

1 

2.5 

3.1 

12.6 

32 

64.6 

106.6 

185 

Conclusion  and  Future  Work:  In  summary,  we  achieved  significant  increase 
in  the  robustness  to  pose  errors  compared  to  [1]  by  considering  all  images  simul¬ 
taneously,  instead  of  using  subsets.  Experimentally,  our  method  ensured  opti¬ 
mality  for  small  pose  errors.  C-arm  tracking  is  a  cumbersome  process  and  easing 
on  this  constraint  has  great  practical  utility.  With  our  method,  a  less  accurate 
estimation  of  the  pose  may  suffice.  For  example  in  [4],  starting  from  an  initial 
guess,  pose  was  further  estimated  iteratively  using  the  current  3D  reconstruc¬ 
tion,  yet  the  seed  matching  remained  susceptible  to  pose  errors.  Applying  our 
method  to  [4]  promises  a  clinically  viable  solution  without  using  external  tracker 
or  encoder  on  the  C-arm.  We  are  also  extending  the  method  to  reconstructing 
overlapping  seeds  that  are  occluded  in  one  or  more  X-ray  images  [5,  6]. 
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ABSTRACT 

C-arm  fluoroscopy  is  modelled  as  a  perspective  projection,  the  parameters  of  which  are  estimated  through  a  cali¬ 
bration  procedure.  It  has  been  universally  accepted  that  precise  intra-procedural  calibration  is  a  prerequisite  for 
accurate  quantitative  C-arm  fluoroscopy  guidance.  Calibration,  however,  significantly  adds  to  system  complexity, 
which  is  a  major  impediment  to  clinical  practice.  We  challenge  the  status  quo  by  questioning  the  assumption  that 
precise  intra-procedural  C-arm  calibration  is  really  necessary.  Using  our  theoretical  framework,  we  derive  upper 
bounds  on  the  effect  of  mis-calibration  on  various  algorithms  like  C-arm  tracking,  3D  reconstruction  and  surgical 
guidance  in  virtual  fluoroscopy  -  some  of  the  most  common  techniques  in  intra-operative  fluoroscopic  guidance. 
To  derive  bounds  as  a  function  of  mis-calibration,  we  model  the  error  using  an  affine  transform.  This  is  fairly 
intuitive,  since  small  amounts  of  mis-calibration  result  in  predictably  linear  transformation  of  the  reconstruction 
space.  Experiments  indicate  the  validity  of  this  approximation  even  for  50  mm  mis-calibrations. 

The  problem  is  twofold:  (a)  C-arm  intrinsic  calibration;  and  (b)  C-arm  distortion  correction.  Using  our 
theoretical  and  experimental  analysis  on  mis-calibrated  C-arms,  we  propose  a  hypothesis  that  indicates  that 
intrinsic  calibration  might  not  be  required  for  a  family  for  surgical  procedures.  The  framework  also  makes 
suggestions  to  the  current  work  flow  so  as  to  further  minimize  the  effect  of  any  inherent  mis-calibration.  To 
address  the  problem  of  pose  dependant  distortion  correction,  we  propose  the  use  of  a  framework  that  can 
statistically  study  the  maximum  variation  in  distortion  near  a  certain  pose  and  then  intra-operatively  use  an 
average  correction.  These  theoretical  derivations  and  experimental  results  make  a  strong  case  for  the  use  of  mis- 
calibrated  C-arms,  obviating  the  cumbersome  intra-operative  calibration  process,  potentially  boosting  clinical 
applicability  of  quantitative  fluoroscopy  in  many  procedures.  We  also  prove  this  hypothesis  experimentally. 

Keywords:  C-arm  Calibration,  3D  fluoroscopic  guidance. 

1.  INTRODUCTION 

C-arm  fluoroscopy  is  ubiquitous  in  general  surgery,  due  to  its  real-time  nature,  versatility,  and  low  cost.  At 
the  same  time,  quantitative  fluoroscopy  has  not  found  a  large  scale  clinical  acceptance ,  because  of  some  inherent 
technical  difficulties.  It  needs  to  solve  four  major  problems:  (1)  C-arm  image  distortion;  (2)  Calibration  of  model 
parameters;  (3)  Pose  recovery  or  tracking  when  multiple  images  are  taken;  and  (4)  Registration  to  other  imaging 
modalities.  Quantitative  fluoroscopy  could  enable  a  significant  improvement  in  the  current  clinical  practice,  the 
prominent  works  tackling  some  of  the  above  problems.1,2 

If  is  known  that  both  image  distortion3  and  intrinsic  calibration4, 5  may  vary  significantly  with  pose.  The  in¬ 
trinsic  parameters  themselves  can  be  classified  into  two  broad  categories  -  pixel  size  and  location  of  X-ray  source 
(focal  spot)  with  respect  to  the  image  intensifier.  The  pixel  size  of  a  C-arm  is  usually  constant  throughout  the 
life  of  the  C-arm,  depending  on  the  actual  electronic  hardware.  The  location  of  the  X-ray  source  on  the  other 
hand  varies  from  pose  to  pose  due  to  the  flexibility  in  the  mechanical  design  and  the  large  weight  of  the  X-ray 
source. 
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Image  distortion  usually  has  a  consequential  contribution  to  reconstruction  error  and  needs  to  be  compen¬ 
sated.  Thus  the  additional  cost  of  a  full  online  calibration  is  not  substantial.  Recently  developed  advanced 
intensifier  tubes  allow  for  lesser  distortion,  while  modern  flat  panel  detectors  obviate  distortion  correction  alto¬ 
gether.  This  fact  brings  up  the  question  whether  we  need  to  calibrate  the  C-arm  fully  at  each  pose.  The  question 
also  leads  to  the  broader  issue,  that  even  if  it  is  not  pose  dependent,  how  accurate  does  calibration  need  to  be. 
In  spite  of  the  importance  of  calibration  in  C-arm  fluoroscopy,  as  far  as  the  authors  are  aware,  there  has  been 
no  prior  work  that  conducts  this  kind  of  analysis.  The  vision  community  has  a  similar  problem6, 7  when  cameras 
are  used  for  visual  serving  of  robots.  In  essence,  this  paper  addresses  the  issue  of  sensitivity  of  3D  fluoroscopy 
to  mis-calibration. 


Figure  1.  The  C-arm  intrinsic  calibration  is  pose  dependant.  Hence  conventional  methods  for  fluoroscopic  guidance  need 
to  calibrate  the  C-arm  at  every  pose,  a  significant  liability  for  intra-operative  navigation.  Traditionally,  a  calibration 
phantom  is  attached  to  the  image  intensifier,  an  X-ray  image  is  taken  and  sent  to  a  calibration  module.  The  phantom  is 
then  removed,  making  the  C-arm  ready  for  use  with  the  patient.  This  is  repeated  for  every  desired  pose. 

In  quantitative  C-arm  fluoroscopy,  we  typically  need  to  measure  the  spatial  transformation  between  two  ob¬ 
jects,  such  as  a  vertebra  and  a  bone  drill,  as  compared  to  the  transformation  between  an  object  and  the  C-arm 
itself.  Thus  the  central  intuition  of  this  paper  is  that  while  an  incorrect  calibration  gives  erroneous  estimates  for 
the  absolute  transformations,  nevertheless  it  still  provides  acceptable  relative  estimates.  The  consequence  of  this 
conjecture  is  potentially  far  reaching,  as  it  can  turn  fluoroscopy  to  an  affordable  quantitative  measurement  tool 
in  a  large  family  of  procedures.  It  should  however  be  noted  that  we  do  not  claim  that  calibration  would  always 
be  unnecessary,  since  there  are  some  applications  that  require  accurate  absolute  estimates,  as  can  be  observed 
in  cone  beam  computed  tomography.5  The  decision  should  always  be  made  case  by  case,  experimentally. 


The  driving  application  for  the  proposed  hypothesis  is  virtual  fluoroscopy  (VF),  where  synthetically  cre¬ 
ated  X-ray  images  in  conjunction  with  image-guidance  aids  the  surgeon  in  accurately  targeting  the  anatomy.8 
Furthermore,  it  significantly  decreases  the  amount  of  radiation  dose  imparted  to  the  surgeon,  which  has  been 
typically  very  high  in  these  family  of  procedures.  In  this  paper,  we  build  a  mathematical  framework  to  formally 
address  this  issue  and  lend  credit  to  the  intuition  that  a  loose  estimate  of  the  C-arm  parameters  might  suffice. 
Furthermore,  our  theoretical  derivations  also  provide  us  with  minor  modifications  to  the  current  VF  procedures, 


so  as  to  minimize  the  effect  of  any  inherent  mis-calibration.  In  essence,  we  prove  in  theory  and  demonstrate 
experimentally  that  VF  is  feasible  with  an  un-calibrated  C-arm. 

2.  MATHEMATICAL  FRAMEWORK 

C-arm  Imaging:  Geometric  aspects  of  fluoroscopic  imaging  can  be  modelled  as  a  perspective  transformation 
with  five  parameters  -  focal  length,  image  origin  and  pixel  size  (Figure  2).  The  transformation  formula  is  given 
in  equation  (1),  (2);  where  Fo,Fx,Fj  are  the  coordinate  frames  of  an  object,  X-ray  source  and  the  image, 
respectively;  P  is  a  3D  point;  Po  are  the  homogenous  coordinates  of  P  in  Fo;  p  is  the  projection  of  P  on  the 
image  plane;  pi  are  the  homogenous  coordinates  (in  pixels)  of  p  in  frame  Fj;  xFo  is  the  4x4  rigid  transformation 
matrix  that  transforms  a  point  in  Fq  to  Fx  (also  called  pose);  1  Fx  is  the  3x4  perspective  projection  matrix;  f 
is  the  focal  length;  O  =  (ox,  oy)  is  the  projection  of  the  X-ray  source  on  the  image  plane  (later  referred  to  as  the 
origin);  sx  and  sy  are  the  pixel  sizes  along  the  X  and  Y  axis  of  the  image. 


Pi  =  1  Fx  X  Ff  Pf  (1) 


Figure  2.  Projective  geometry  and  notations  for  fluoroscopic  imaging. 

There  are  a  total  of  five  independent  parameters  that  need  to  be  evaluated  by  the  calibration  procedure  - 
the  pixel  sizes  (two)  and  the  focal  spot  (three).  Note  that  owing  to  the  mathematics  of  projective  imaging, 
these  intrinsic  parameters  can  only  be  estimated  up  to  an  arbitrary  scale  (which  determines  the  physical  size  of 
the  C-arm).  Hence  instead  of  the  pixel  sizes,  the  aspect  ratio  (ratio  of  pixel  sizes)  is  typically  used.  Since  the 
aspect  ratio  remains  unchanged  throughout  the  life  of  the  C-arm,  online-calibration  essentially  reduces  just  to 
the  location  of  the  focal  spot.  The  proposed  theoretical  framework  can  be  used  to  study  sensitivity  due  to  any 
of  the  five  parameters,  nevertheless,  we  limit  ourselves  only  to  that  of  the  focal  spot. 

2.1.  Model  for  Reconstruction  Space  Transformation 

As  illustrated  in  Figure  8,  let  A  &  B  (with  reference  frames  Fa  &  Fb)  be  the  two  objects  being  imaged.  The 
assumptions  are:  (i)  1  Fa/  Fb  can  be  computed  from  the  images;  (ii)  A  &  B  are  not  large  in  comparison  to  the 
focal  length;  (iii)  Fa  and  FB  are  close  by ;  and  (iv)  the  quantity  of  interest  is  AFB  =  (7Fa)_1  1  Fb.  Let  /i  be 
the  true  focal  spot  and  =  (/i  +  D)  be  the  mis-calibrated  estimate.  We  claim  that  even  though  the  absolute 
locations  of  the  objects  are  off,  their  relative  transformation  might  still  be  accurate. 


['FaV i)]-1  7Fb(/i) 


['Ci  (/-»)]  1  7Fb(/2) 


(3) 


Wrong  X-ray 
source 


Figure  3.  Mis-calibration  shifts  all  reconstructed  objects  under  an  affine  transformation. 

A  transformation  is  needed  that  can  take  the  absolute  location  of  an  object  reconstructed  with  calibration 
/i,  and  compute  its  corresponding  location  with  calibration  fa.  We  claim  that  the  simplest  transformation  will 
be  a  linear  affine  model  T.  The  intuition  derives  from  the  observation  that  the  image  plane  is  the  same  in  both 
reconstruction  spaces.  Thus  if  Pi  (not  in  homogenous  coordinates)  projects  to  a  point  p  on  the  image,  then 
it  is  constrained  to  be  on  line  L\  in  the  /i-space  and  on  L 2  in  /2-space.  Thus  we  seek  a  continuous  invertible 
transformation  that  projects  L\  to  £2-  By  incorporating  the  above  constraints,  T  can  be  evaluated  to  be, 


P2  =  T  -Pi 


'1  0  Dx/flz 

0  1  Dy/flz 

0  0  1  +  (Dz/flz) 


• Pi  =  Pi  +  (d-Z/flz)D 


(4) 


where  with  respect  to  (wrt)  F/,  D  =  (Dx,Dy,Dz)]  d  =  ||P||2;  D  =  D/d ;  /1  =  (flx,  fly,  flz);  and  Pi  = 
(X,  Y,Z).  Each  point  is  effectively  translated  in  direction  D  by  an  amount  proportional  to  its  distance  from  the 
image.  Experiments  measuring  the  correctness  of  this  affine  model  are  available  in  Section  3.  Thus  to  study 
sensitivity,  it  is  sufficient  to  study  the  properties  of  T. 


Changes  in  Length  and  Scale  T  preserves  the  scale  along  the  x,  ?/-axes,  but  scales  the  space  along  the  z- 
axis.  Let  Pi (Xi,  Yi,  Zi)  &  P2(X2,  Y2,  Z2)  be  any  two  points  (not  necessarily  close  to  each  other)  in  the  /i-space 
at  a  distance  of  l±.  T  maps  them  to  points  Q 1  &  Q 2  in  the  /2-space  at  a  distance  of  I2  (Figure  4  (a)).  It  can  be 
shown  that 


l22  =  li  +  i-(Zi  -  z2)2  +  2 (Zl  Zz)  D- (Pi-  P2)  (5) 

Jlz  Dz 

Wh  -/ill  <  P\Z1-Z2\  (6) 

Jlz 

=>  h[l-P]  <  h  <  h[l  +  P]  (7) 

Jlz  Jlz 

It  directly  follows  from  Equation  (6)  that  T  does  not  alter  the  length  significantly.  As  an  example,  a  10  mm 
calibration  error  would  affect  the  length  of  a  30  mm  thoracic  pedicle  screw  at  an  angle  of  45°  by  less  than  0.2  mm 
(focal  length  ~  1  m),  which  is  significantly  less  than  the  error  from  other  sources.  Thus  Fa,Fb  will  not  change 
their  relative  translation  by  a  factor  more  than  that  specified  by  Equation  (6). 


Figure  4.  (a)  Mis-calibration  rotates  and  scales  a  straight  line  segment,  (b)  Pose  dependent  calibration  might  be 
successfully  approximated  by  using  the  mean  value. 


Changes  in  Absolute  Orientation  A  change  in  orientation  results  from  the  object  having  a  depth  (Figure 
4  (a)).  It  can  be  shown  geometrically  that  the  orientation  error  is  maximal  when  the  vector  P1P2  is  roughly 
orthogonal  to  D  and  is  purely  in  the  vertical  plane.  The  amount  ( 0 )  and  the  axis  (k)  of  rotation,  through  a  series 
of  computations  can  be  shown  to  be  as  in  Equation  (8).  The  bound  on  the  rotation  error  is  dependent  only  on 
origin  mis-calibration  and  not  on  that  in  focal  length  .  More  importantly  it  is  independent  of  the  height /depth  of 
the  object  (as  far  as  it  is  non-planar)  and  its  distance  from  the  image  plane.  Thus  Fa,  Fb  in  Figure  8  will  observe 
the  same  absolute  rotation,  in  effect  not  experiencing  any  relative  rotation.  Experimental  results  corroborating 
this  claim  are  available  in  Section  3. 
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Error  in  Reconstruction  of  Point  Features  In  many  applications  (particularly  in  brachytherapy ) ,  C-arms 
are  used  to  reconstruct  3D  point  objects.  This  is  done  by  obtaining  multiple  images  at  varying  orientations  and 
then  using  triangulation  to  obtain  the  desired  intersection.  In  ideal  circumstances,  all  the  lines  would  intersect 
at  a  unique  point.  In  practice  however,  calibration  (and  other)  errors  lead  to  non-intersecting  lines.  We  will 
attempt  to  bound  the  error  in  this  symbolic  reconstruction  of  the  point.  Let  point  P  be  imaged  from  N  different 
poses  and  reconstructed  in  a  tracker  frame  Fa,  which  is  stationary  wrt  P.  Let  the  ith  pose  have  a  focal  spot 
error  (in  frame  Fa)  of  Di.  Without  errors,  each  reconstructed  line  (li)  would  pass  through  P.  It  can  be  shown 
that  due  to  the  calibration  error  ZA,  the  new  line  passes  through  a  new  point  P'A  and  undergoes  a  rotation  <fr. 
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Pa  0  0 


(PA-Pj), 

fiz  1 


( Ji  '  Pi)  sinOi 
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(9) 


where  0i  is  the  angle  that  li  makes  with  the  z-axis  of  Fa-  The  rotation  is  fairly  small  and  can  be  ignored.  Thus 
Pa  is  at  a  distance  of  (Pa  •  Di) sinOi / fiz  from  li.  If  Q  is  the  symbolic  intersection  of  all  li  s,  then  it  can  be  shown 
that  Q  is  no  further  away  than  ( d™ax  sinOmax )  \ \ Pa  \ |  away  from  any  of  the  lines.  Moreover,  the  reconstruction 
error  (RE)  can  also  be  shown  to  be  bounded  by 


RE  =  \\(Q  -  PA)\\  <  ||C4|| 

Jz 


(10) 


where  dmax  is  the  maximum  amount  of  mis-calibration  and  fz  is  the  minimum  focal  length.  Thus  a  10  mm 
focal  length  error  causes  an  error  less  than  0.5  mm  for  a  point  at  a  distance  of  35  mm.  Note  that  this  is  the 
worst  case  error  analysis  and  in  practice  the  dot  product  in  Equation  (9)  mutually  cancels  positive  and  negative 
errors,  leading  to  extremely  low  reconstruction  errors  (Section  3). 


The  Optimal  Choice  for  Calibration  Parameters  Since  the  focal  spot  is  pose  dependant,  and  the  previ¬ 
ous  results  suggest  robustness  to  mis-calibration,  choosing  a  constant  calibration  for  quantitative  reconstruction 
might  be  viable.  In  the  scenario  that  the  focal  spot  might  vary  as  much  as  10  mm  from  one  pose  to  another, 
“ what  constant  calibration  should  be  chosen  to  minimize  error”! 

Let  us  assume  that  we  are  imaging  a  point  P  from  N  different  poses  (Figure  2  (c)).  Wrt  frame  Fj,  let  the  ith 
pose  have  the  focal  spot  at  fa  =  (fax,  fay ,  faz )  and  the  point  be  at  location  Pi  =  (JQ,  Zi).  Note  that  we  assume: 
(a)  variations  in  each  of  fax,  fay,  fiz,Xi,Yi,  Zi  &  pose  are  independent;  (b)  Pi  s  are  close  to  the  iso-center,  i.e. 
variations  in  Xi,Yi,Zi  are  not  high.  We  choose  a  constant  value  of  F  =  (Fx,Fy,Fz)  for  the  focal  spot,  which 
will  displace  the  point  Pi  to  Qi  =  T(fa,F)  • Pi .  The  aim  is  to  choose  an  F  which  minimizes  the  net  variation  of 
A Qi  —  Qi  —  hq.  Through  a  series  of  computations,  it  can  be  shown  that 


A  Qi  =  (Pi  -  up) 


MQ  =  Mp  +  —(F  -  [if) 
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where  /jlq,  /ip,  /jlz,  fjbfz,jlf  are  the  mean  values  of  Qj,  Pj,  Zj,  fjz,  fa;  Zj  =  fiz  +  AZj  and  likewise  for  fjz,fa, 
where  j  =  1 . . .  N.  In  the  above  calculations,  the  second  order  terms  either  summed  to  0  due  to  the  independence 
of  the  variables  or  were  too  small  in  comparison.  Our  choice  of  F  should  be  the  one  that  minimizes  the 
variance(AQ)  =  var(AQx)  +  var(AQy)  +  var(AQz).  It  should  be  noted  that  Fz  scales  the  whole  space,  i.e. 
a  lower  value  will  decrease  the  variance,  implying  that  the  choice  of  Fz  =  0  forces  var(Qz)  =  0  by  forcing  all 
Q[s  to  lie  on  a  plane.  Thus  var(Qz )  does  not  provide  sufficient  constraints  for  Fz.  We  will  first  obtain  FX1  Fy 
by  minimizing  the  variance  along  x,y- axes  (since  there  is  no  scaling  in  these  directions),  and  then  will  compute 
Fz.  Notice  that  the  first  term  in  Equation  (12)  is  due  to  the  relative  movement  in  P,  while  the  rest  is  due  to  an 
error  in  the  calibration.  Since  we  are  interested  only  in  the  variance  due  to  mis-calibration,  we  will  ignore  the 
variations  in  P.  Minimizing  var(AQ)  and  enforcing  independence  of  fax,fay  &  faz  gives 


F  = 


AfazAfj 

SfA  fl 


h fz  —  [  llfx  5  M/y?  0  ] 
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As  expected,  Fz  =  0  from  above.  To  compute  Fz,  we  need  to  impose  length  preserving  constraints.  Thus 
if  we  measure  a  line  segment  of  length  l  in  each  image,  use  Equation  (6)  to  derive  the  net  length  error,  the 
minimization  implies 


Fz 
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(14) 


Thus  F  =  ft f  (the  mean),  which  is  fairly  intuitive  and  probably  in  common  practice.  Likewise,  this  particular 
choice  of  Fx,  Fy  is  also  a  length  preserving  constraint,  i.e.  it  minimizes  the  error  in  lengths  of  line  segments. 
Calibration  error  in  A  Qi  now  reduces  to  —-jjj-Afa,  which  has  a  stable  mean  and  low  variance.  Equation  (15) 

gives  a  bound  on  the  error  when  the  assumed  value  of  F  is  away  from  the  mean  jlf  by  a  distance  d.  A  10  mm 
variation  in  the  focal  length  (var  ~  3  mm),  P's  roughly  at  the  iso-center  having  a  depth  variation  of  100  mm 
and  the  assumed  calibration  unusually  away  from  \if  by  50  mm  still  bounds  the  maximum  error  by  0.75  mm. 


Thus  large  and  constant  mis-calibration  in  many  applications,  might  still  provide  sub-millimetric  3D  quantitative 
measurements. 


error  < 


y/d2var(Z)  +  n2zvar{\\f\\) 
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2.2.  Mis-calibration  and  its  effect  on  Virtual  Fluoroscopy 

In  traditional  C-arm  fluoroscopy,  surgeons  often  use  C-arms  whenever  they  want  to  know  the  location  of  a  tool 
inside  the  patient.  In  virtual  fluoroscopy  (VF),  X-ray  images  are  created  synthetically  with  the  position  of  the 
tool  (tracked  in  3D)  superimposed  on  the  2D  images.  This  not  only  adds  to  the  targeting  accuracy,  but  also 
reduces  the  radiation  dose  both  to  the  surgeon  and  the  patient.  Thus  VF  allows  the  tool’s  position  to  be  known 
without  taking  a  new  image  each  time.8 


(a) 


(b) 


Figure  5.  (a)  The  overall  dataflow  for  a  virtual  fluoroscopy  system.  The  current  3D  tool  location  is  registered  to  a 
previously  taken  C-arm  image  of  the  patient,  and  then  rendered  artificially  on  the  image.  Two  such  images  can  provide 
real-time  3D  guidance,  without  the  need  to  take  more  images,  (b)  An  example  screen-shot  of  our  implemented  virtual 
fluoroscopy  system.  The  ’x’  denotes  the  system’s  projected  tool-tip  position,  overlaid  on  top  of  the  corresponding  dewarped 
image.  For  reference,  the  true  location  of  the  projected  tool  tip  is  also  visible  in  the  X-ray  image.  The  FTRAC  fiducial 
used  to  register  the  tracker  to  the  C-arm  is  also  visible. 


As  illustrated  in  Figure  5,  the  tool  is  tracked  realtime  in  3D  using  an  auxiliary  tracking  system.  The  tracker 
is  registered  to  the  C-arm  using  some  sort  of  a  prefabricated  precision  phantom,  visible  both  to  the  tracker  and 
in  the  X-ray  image.  This  allows  for  the  projection  of  the  current  3D  tool  location  on  a  previously  taken  X-ray 
image  of  the  anatomy,  essentially  creating  a  ’virtual’  X-ray  image.  A  prerequisite  for  projecting  the  current  tool 
location  on  the  image  is  to  know  the  location  of  the  X-ray  image  wrt  the  X-ray  source  accurately  (otherwise 
known  as  C-arm  calibration). 

Traditionally,  the  C-arm  is  extensively  calibrated  with  each  change  in  pose  before  taking  an  image,  which 
is  a  cumbersome  procedure.  A  typical  such  phantom  is  illustrated  in  Figure  1.  This  naturally  leads  us  to 
ask  ’is  it  actually  necessary  to  calibrate  the  C-arm  for  every  change  in  pose?’.  More  importantly,  how  precise 
does  the  calibration  need  to  be  for  a  desired  accuracy  in  a  surgical  procedure.  In  what  follows,  using  the  mathe¬ 
matical  framework  developed  so  far,  we  study  the  effect  of  this  mis-calibration  on  the  accuracy  of  tool  projection. 

To  develop  the  framework,  we  constrain  the  tool  only  to  a  point,  which  is  reasonable  since  most  tools  are  a 
collection  of  3D  points  attached  rigidly  to  the  tracker  on  the  tool  (for  ex.  an  insertion  needle).  We  start  with 
adding  a  calibration  error  of  D  =  (Dx,  Dy,  Dz)  in  the  source  position  wrt  the  image.  The  results  in  Equation 


(4),  (8)  indicate  that  the  primary  errors  are  in  the  plane  containing  the  true  X-ray  source,  incorrect  X-ray  source 
and  orthogonal  to  the  image.  We  will  refer  to  this  plane  as  the  error-plane  and  is  illustrated  in  Figure  6.  Thus 
this  plane  can  be  assumed  to  be  the  base  plane  with  the  mis-calibration  D  =  +  D 0,  Dz).  For  simplicity 

of  calculation,  we  will  drop  the  y-coefhcient  and  assume  D  =  L/D%  +  D^Dz). 


Projection  Error 


Figure  6.  With  the  true  calibration  /i,  the  3D  point  Pi  wrt  F\  projects  to  pi  in  the  image.  Due  to  the  incorrect 
estimate  /2,  Fi  is  mis-estimated  at  Fh,  moving  the  point  from  Pi  to  P2,  eventually  projecting  at  P2  in  the  image.  P3  is  a 
convenient  3D  point,  which  when  projected  from  / 2  projects  on  pi.  \\pi  —  P2W  represents  the  error  in  the  system  due  to 
mis-calibration.  Note  that  the  error  is  assumed  in  the  X-Z  plane  without  any  loss  in  generality. 

Let  Fi  be  the  true  frame  coordinate  of  the  Tracker  to  C-arm  registration  phantom,  and  F 2  be  the  frame  of 
the  phantom  as  computed  using  an  incorrect  calibration.  Though  illustrated  such,  Fi ,  F \  are  not  necessarily  in 
the  error  plane.  The  relationship  between  F2  and  F\  is  as  evaluated  in  Equation  (4),  (8),  involving  a  translation 
in  the  direction  of  D  and  a  rotation  along  the  Y-axis.  Let  Pi  =  (a,  b )  be  the  true  location  of  the  tool  tip  wrt  the 
phantom  in  frame  F/.  Note  that  this  value  (when  expressed  in  frame  Fi)  is  independent  of  C-arm  calibration 
and  is  computed  using  the  tool  to  tracker-base  transformation  and  tracker-base  to  phantom  registration.  Let  P2 
(expressed  in  frame  Fj)  be  the  new  location  of  the  tool  tip  wrt  the  phantom  when  the  frame  of  the  phantom 
has  been  reconstructed  using  a  mis-calibrated  C-arm  (F2).  Let  pi  be  the  true  projection  of  the  tool  tip  from  Pi , 
and  P2  the  mis-calibrated  projection  of  the  tip  from  P2.  Note  that  pi,P2  will  not  be  the  same  in  general,  the 
difference  representing  the  projection  error  due  to  mis-calibration  in  a  VF  system.  We  will  attempt  to  compute 
an  upper  bound  on  this  error. 


To  compute  ||pi  — P2 1| ,  we  introduce  another  3D  point  P3,  such  the  P3  when  projected  from  the  mis-calibrated 
X-ray  source,  projects  on  the  2D  point  pi.  Thus  points  P2,  P3  when  projected  using  the  mis-calibrated  parameters, 
can  be  used  to  estimate  ||pi  —  _p2 1 1  •  One  example  point  for  P3,  computed  using  Equation  (4),  is  as  satisfied  by 
the  equation  below,  where  Zi  is  the  Z-coordinate  of  point  Pi  wrt  Fj. 

P3~Pi=(j)D  (16) 

Note  that  P2  wrt  frame  F2  is  same  as  Pi  wrt  frame  F\.  Let  this  vector  be  v.  If  6  as  computed  by  Equation 
(8)  be  the  angle  of  rotation  (along  the  new  Y-axis),  then  v  is  rotated  by  6  along  the  Y-axis  of  the  C-arm.  Thus 
P2  can  be  computed  from  Pi  by  applying  the  translation  and  rotation  of  from  Fi  to  F2  as  given  by  Equations(4), 
(8).  If  Zf  is  the  Z-coordinate  of  point  Fi  wrt  F/,  then 


P2  —  P\  =  Translation(y )  +  Rotation(v ) 

=  [(-^-)D\  +  [(b, —a)sinO  —  (a,b)(l  —  cosO)\ 
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Note  that  we  have  not  written  down  the  components  of  Pi,  P2,  P3  along  the  Y-axis.  This  is  because,  the  affine 
transformation  model  does  not  make  any  transformations  (rotation  or  translation)  along  this  axis,  resulting  in 
no  change  from  Pi  to  P2,  P3.  From  Equations  (16)  &  (17),  we  get 


-bD  +  (b,-a)^D2x  +  D2y 

P2-P3  =  7 

-b^Dl  +  DlM  +  ( 6,— a)  ^£>g +  £>2 
=  7  (18) 

aJof^D2  +  bDz 

=  (0’“ - - 7 - } 

=  (0,/3) 

Note  that  the  x-coordinate  of  the  vector  P2  —  P3  is  zero.  We  will  use  this  special  structure  to  evaluate 
||_P2  —  Pi||-  For  this  computation,  we  can  assume  without  loss  of  generality  that  Fj  is  located  at  the  origin  of  the 
mis-calibrated  image  plane.  If  wrt  F/,  P3  =  (x,y,z)  (y-coordinate  included)  then  pi  =  (/'//'  —  z)(x,y),  where 
f  is  the  mis-calibrated  focal  length.  Since  P2  =  P3  +  (0,0,/?),  it  can  be  shown  that 
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where  m  is  the  magnification  factor  of  the  3D  point  and  <f>  is  the  angle  that  the  3D  point  makes  with  the 
projection  axis  of  the  C-arm  (the  line  joining  the  X-ray  source  and  the  origin,  usually  very  close  to  the  center  of 
the  image).  Thus  the  error  due  to  mis-calibration  in  a  virtual  fluoroscopy  system  is 
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From  Equation  (20)  it  can  be  noticed  that  the  error  in  a  VF  system  can  be  significantly  reduced  by:  (a) 
bringing  the  tracked  tool  closer  to  the  center  of  the  image;  and/or  (b)  using  a  registration  phantom  that  can 
work  near  the  tool  workspace.  To  get  a  sense  of  the  numbers,  if  the  tool  is  somewhere  in  the  middle  of  the 
C-arm  workspace  (magnification  =  2)  &  projects  at  the  corner  of  a  9”  image  intensifier  (tancj)  ~  0.12),  then 
a  mis-calibration  as  high  as  even  50  mm  will  not  induce  an  error  greater  than  1  mm  in  a  VF  system,  if  the 
registration  phantom  is  used  inside  a  100  mm  region  of  the  tool  workspace.  Alternately,  even  if  the  registration 
phantom  needs  to  be  kept  far  away  (at  a  distance  of,  say,  250  mm),  then  keeping  the  projection  of  the  tool 
in  a  40  mm  diameter  region  near  the  image  center,  will  still  produce  an  error  less  than  1  mm  with  a  100  mm 
mis-calibration.  Note  that  a  1  mm  error  in  the  image  at  a  magnification  of  2,  indicates  only  a  0.5  mm  3D  error 
(when  multiple  images  are  used,  as  is  customary).  Also  note  that  in  most  cases,  typical  mis-calibrations  are 
predominantly  in  the  imaging  direction,  in  which  case,  it  suffices  if  the  registration  phantom  can  be  used  at  the 
same  ’depth’  as  the  tool  workspace.  In  such  cases,  the  dot  product  in  Equation  (20)  cancels  out,  ensuring  that 
no  error  is  added  in  the  tool  guidance  of  the  VF  system. 

2.3.  Distortion  Correction 

C-arm  images  suffer  from  image  distortion  (as  high  as  5  mm),  owing  to  the  curved  nature  of  detector  (radial 
distortion)  and  the  earth’s  magnetic  field3  (s-shaped  distortion).  In  order  to  achieve  quantitative  surgical  navi¬ 
gation,  the  distortion  needs  to  be  corrected  for.  This  is  cumbersome,  more  so,  since  the  warp  is  pose  dependent. 
Traditional  distortion  correction  techniques  use  a  calibration  phantom1  4  and  characterize  the  distortion  by  fit¬ 
ting  a  high-order  polynomial  between  the  distorted  image  and  the  expected  image.  Though  easy  to  use,  they 
are  cumbersome  for  intra-operative  correction  due  to:  1)  the  interference  between  the  patient  anatomy  and  the 
phantom;  and  2)  increased  radiation  dose  &  surgical  time  when  two  images  are  used,  one  with  and  the  other 
without  the  phantom.  Alternately,  the  distortion  parameters  at  a  particular  pose  could  be  interpolated  from  the 
multiple  neighboring  poses  taken  pre-operatively,3  but  require  accurate  c-arm  tracking.  The  intra-operative  use 
of  the  same  c-arm  pose,  as  was  used  pre-operatively  for  distortion  correction,  has  also  been  proposed,2  though  it 
requires  the  c-arm  pose  to  be  repeatable  and  hence  could  become  time-consuming  when  collecting  many  images. 
All  in  all,  online  distortion  correction  has  been  a  cumbersome  problem,  without  a  reliable/robust  solution. 

In  line  with  the  general  theme  of  the  paper,  we  observe  that  image  distortion  does  not  radically  change 
between  two  neighboring  poses,  i.e.  small  movements  if  the  C-arm  result  in  small  variations  in  image  distortion. 
This  suggests  that  it  might  suffice  to  use  a  constant  distortion  correction  in  a  small  region  of  interest  (ROI). 
Unfortunately,  the  definition  of  a  small  ROI  or  neighboring  poses  is  not  clear  and  tools  that  can  crisply  charac¬ 
terize  these  ideas  have  not  yet  been  developed.  In  what  follows,  we  propose  some  simple  yet  powerful  tools  to 
further  reinforce  these  ideas  from  a  quantitative  perspective. 


Statistical  Characterization  of  C-arm  Distortion:  C-arm  distortion  has  been  traditionally  modelled  using 
a  high-order  polynomial.  This  polynomial  is  computed  by  using  a  phantom  grid  with  radio-opaque  beads  in  a 
known  pattern,  where  the  observed  distorted  image  is  compared  to  the  expected  image.  Among  the  many 
basis  polynomials,  the  authors  have  found  the  Bernstein  basis  of  order  4/5  as  sufficiently  robust  &  fast,  as  also 
observed  in  the  literature.2  Once  the  polynomials  are  computed,  they  can  be  used  to  estimate  the  amount  of 
sub-pixel  distortion  at  any  given  pixel  in  the  image.  This  computation  can  be  done  independently  for  multiple 
pre-operative  images.  To  study  precisely  the  extent  of  statistical  variations  in  distortion  among  images  of  a 
certain  ROI,  a  principal  component  analysis  (PC A)  can  be  conducted  on  all  the  distortion  maps.  To  do  so,  the 
distortion  map  for  each  image,  discretised  at  each  pixel,  is  organized  as  a  single  vector  in  a  very  high  dimensional 
space.  Thus  each  distortion  map  is  a  single  point  in  this  high-dimensional  space.  To  quantitatively  study  the 
amounts  of  variation,  it  shall  suffice  to  study  the  statistics  of  this  collection  of  points.  The  PC  A  analysis  can  be 
of  tremendous  aid  here,  and  can  provide  a  linear  basis  for  the  subspace  of  these  points,  as  given  below. 
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where  D  is  the  distortion  at  any  given  C-arm  pose,  M  is  the  mean  distortion  map  of  the  whole  space  and 
Di  is  the  ith  strongest  (principal)  component  (mode).  Furthermore,  the  eigen- values  (computed  during  PC  A) 
associated  with  each  principal  mode  indicate  the  variance  along  each  associated  direction.  This  allows  for  a 
statistical  analysis  of  the  pose-dependant  C-arm  distortion.  It  has  been  observed  experimentally  that  more  than 
98%  of  the  distortion  can  be  modelled  using  just  the  first  three  principal  modes,  indicating  that  though  complex, 
pose  dependant  C-arm  distortion  lies  predominantly  inside  a  3-dimensional  linear  space.9 

This  experimentally  determined  numerical  observation  (for  a  particular  C-arm)  that,  a  small  movement  in  the 
C-arm  pose  will  result  in  only  a  small  variation  in  image  distortion,  can  be  used  to  motivate  two  simple  algorithms 
for  distortion  correction:  (a)  use  the  mean  computed  from  Equation  (21)  for  correcting  the  distortion  within  a 
certain  acceptable  region;  (2)  if  the  mean  computation  is  unfeasible,  use  a  center  image  that  is  geographically 
close  to  the  mean  image.  Intuitively,  if  an  intra-operative  C-arm  image  is  desired  near  a  certain  pre-calibrated 
pose,  then  compute  a  mean  distortion  near  that  particular  pose.  Furthermore,  from  a  quantitative  3D  guidance 
perspective,  the  use  of  the  pre-operative  distortion  map  itself  might  suffice  in  some  cases. 

3.  PHANTOM  EXPERIMENTS  AND  RESULTS 

Validity  of  the  Model:  Equations  (4)  &  (8)  give  the  translation  and  rotation  transformations  as  predicted  by 
the  affine  model,  the  accuracy  of  which  would  furnish  the  validity  of  the  model.  We  used  the  FTRAC  fiducial 
(Figure  7),  a  small  image-based  fluoroscope  tracking  fiducial,  which  (given  the  calibration)  can  track  a  C-arm 
with  an  accuracy  of  0.5  mm  in  translation  and  0.65°  in  rotation.10  The  fiducial  was  imaged  using  a  Philips 
Integris  V3000  fluoroscope  and  the  true  calibration  read  off  the  machine  display.  The  images  were  not  corrected 
for  distortion.  The  pose  of  the  fiducial  (wrt  to  Fj)  was  first  evaluated  using  the  correct  calibration,  and  then  with 
the  mis-calibrated  parameters.  The  difference  between  the  pose  change  predicted  by  the  equations  and  the  one 
computed  using  the  non-linear  pose  estimation  software,  is  displayed  in  Figure  8  (a)  as  a  function  of  maximum 
calibration  error.  A  total  of  40,  500  points  (from  total  75  poses)  with  randomly  added  mis-calibration  were  used. 
Even  when  mis-calibration  is  as  high  as  25  mm,  the  model  can  predict  the  rotation-axis  with  an  accuracy  of 
4°,  amount  of  rotation  by  0.1°  and  translation  under  5  mm.  For  extreme  mis-calibrations  the  translation  error 
linearly  increases,  while  rotation  is  still  stable.  Thus  the  model  seems  to  predict  with  an  acceptable  accuracy. 


Figure  7.  An  image  of  the  seed  phantom  attached  to  the  FTRAC  fiducial  (left).  The  seed  phantom  can  replicate 
any  implant  configuration,  using  the  twelve  5  mm  slabs  each  with  over  a  hundred  holes.  A  typical  X-ray  image  of  the 
combination  (right). 


Accuracy  of  C-arm  Tracking:  The  FTRAC  fiducial  was  mounted  on  a  0.02°  accuracy  rotational  turntable, 
while  the  fluoroscope  was  kept  stationary.  The  turntable  was  rotated  by  known  precise  amounts  (ground  truth) 
and  images  were  taken.  The  relative  poses  were  also  computed  using  the  pose  estimation  software.  The  accuracy 
in  the  estimation  of  C-arm  motion  is  given  by  the  difference  between  the  computed  relative  pose  and  the  true 
relative  pose.  The  tracking  accuracy  is  plotted  in  Figure  8  (b)  as  a  function  of  mis-calibration.  Even  a  high 
mis-calibration  of  150  mm  adds  no  additional  error  in  C-arm  motion  estimation,  fixing  the  value  at  0.45  mm  in 


translation  and  0.6°  in  rotation.  An  unusually  high  mis-calibration  of  400  mm  also  only  marginally  decreases 
accuracy.  Thus,  mis-calibration  does  not  increase  the  error  of  C-arm  tracking  . 


Error  in  the  C-arm  Source  Position  (mm) 


Sensitivity  of  Relative  C-Arm  Pose  Estimation  to  Calibration 
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Figure  8.  Note  the  scale  variation  in  x-axis.  (a)  An  affine  transformation  is  able  to  predict  the  movement  of  3D  objects 
due  to  mis-calibration;  (b)  C-arm  tracking  is  insensitive  to  mis-calibration;  3D  Reconstruction  is  insensitive  to  mis- 
calibration  in  (c)  origin  ;  (d)  focal  length  up  to  50  mm,  beyond  which  it  starts  to  linearly  drift  away  from  the  tracking 
fiducial.  Notice  that  the  shape  of  the  implant  (relative  err)  is  barely  altered;  (e)  3D  reconstruction  error  decreases  with 
an  increase  in  images  used. 


3D  Quantitative  Reconstruction  using  Multiple  Images:  In  addition  to  tracking  a  C-arm,  it  is  equally 
important  that  multiple  objects  in  the  field  of  view  (eg.  vertebrae  and  screws)  be  reconstructed  accurately 
relative  to  each  other.  In  order  to  validate  our  hypothesis  that  3D  reconstruction  might  not  be  sensitive  to 
mis-calibration,  we  use  an  accurate  acetol  phantom  (Figure  7)  having  100  dummy  radioactive  seeds,  approxi¬ 
mating  a  brachytherapy  implant  (Figure  7).  The  true  3D  coordinate  of  each  seed  wrt  the  fiducial  is  known  by 
rigid  attachment.  The  C-arm  is  tracked  using  the  FTRAC  fiducial  and  the  3D  seed  coordinates  are  computed 
by  triangulation  (an  algorithm  called  MARSHAL  is  used  to  establish  correspondences).  The  difference  between 
the  computed  and  the  true  seed  location  gives  us  the  3D  reconstruction  error  for  each  seed  (wrt  fiducial).  The 
relative  reconstruction  error  removes  any  consistent  shift  reflecting  any  change  in  shape.  These  errors  are  plotted 
as  a  function  of  mis-calibration  in  Figure  8  (c),  (d).  The  reconstruction  error  is  insensitive  to  mis-calibration 
in  origin  and  focal  length  errors  of  up  to  50  mm.  The  shape  of  the  implant  is  stable  even  for  large  calibration 
errors.  Figure  8  (e)  shows  a  drop  in  reconstruction  error  as  the  number  of  images  increase.  Thus  mis-calibration 
does  not  decrease  reconstruction  accuracy. 

3D  Virtual  Fluoroscopy:  We  implemented  an  experimental  setup  in  order  to  evaluate  the  sensitivity  of  a 
VF  system  to  C-arm  mis-calibration,  the  main  components  of  the  system  being:  (a)  A  GE  9600  C-arm;  (b)  a 
precise  C-arm  distortion  correction  and  calibration  toolkit  (implemented  in  Matlab);  (c)  a  Polaris  tracker;  (d)  a 
calibrated  Polaris  pointer  (calibrated  using  the  NDI  ToolViewer);  (e)  the  FTRAC  fiducial. 


Distortion  correction  and  precise  C-arm  calibration  is  performed  for  the  C-arm  images  using  the  phantom 
shown  in  Figure  1.  The  phantom  is  a  typical  two  plate  contraption  that  can  be  attached  to  the  intensifier.  The 
bottom  plate  with  the  evenly  spaced  metal  BB’s  corrects  for  any  image  distortion.  Each  straight  line  on  the 
upper  plate,  along  with  its  projected  image-line  constraints  the  X-ray  source  to  a  3D  plane.  The  least  square 
intersection  of  these  ten  planes  provides  an  accurate  estimate  of  the  source  location  wrt  the  image  (C-arm  cali¬ 
bration).  In  order  to  display  the  current  location  of  the  tool  in  a  VF  system,  the  tracker  (Polaris)  needs  to  be 
registered  to  the  C-arm.  By  virtue  of  its  design,  the  FTRAC  fiducial  can  register  itself  to  the  C-arm  (in  every 
image  that  it  is  present).  Thus,  the  Polaris  to  C-arm  registration  can  be  achieved  by  registering  the  FTRAC 
fiducial  to  the  Polaris.  To  realize  this,  the  tracked  Polaris  pointer  was  used  to  digitize  4  straight  lines  on  the 
FTRAC  fiducial,  the  exact  location  of  which  was  precisely  known  from  the  design  file.  A  non-linear  optimization 
algorithm  aligns  these  straight  lines,  essentially  yielding  the  FTRAC-Polaris  transformation.  Note  that  this 
registration  is  C-arm  pose  independent. 

At  any  given  pose  of  the  C-arm,  the  calibration  phantom  is  used  to  dewarp  the  image  and  to  compute  the 
true  calibration.  The  phantom  is  then  detached  and  a  new  image  is  taken  with  both  the  pointer-tip  and  the 
FTRAC  fiducial  visible  in  the  image  (attached  steadily  using  clamps) .  The  3D  pointer  location  is  measured  using 
the  Polaris.  The  location  of  the  pointer  tip  can  now  be  projected  in  the  X-ray  image  using  the  transformation 

Pproj  =  TFc  *  C Ff  *  F Fp  •  pFt  •  ptip  (22) 

where  pup  is  the  location  of  the  tool-tip  wrt  the  tool  frame  (Ft),  pFt  is  the  measured  transformation  from 
the  tool  frame  to  the  Polaris  frame  (Fp),  F Fp  is  the  pre-computed  transformation  from  the  Polaris  frame  to  the 
FTRAC  frame  (Fp),  c Ff  is  the  transformation  from  the  FTRAC  to  the  C-arm  frame  (Fc)  computed  for  every 
image,  IVc  is  the  projective  transformation  that  converts  any  3D  point  to  its  projection  in  the  2D  image  (Fj), 
and  Pproj  is  the  (computed)  VF  projected  tool-tip  location  in  the  image.  Note  that  IVc  and  c Fp  will  change 
from  one  C-arm  mis-calibration  to  another.  If  pseg  is  the  true  projection  of  the  tool-tip,  segmented  from  the 
X-ray  image,  then  \\pproj  —  PsegW  provides  the  Euclidean  error  estimate  in  the  VF  system. 

A  total  of  75  such  X-ray  images  were  acquired  by  randomly  varying  the  C-arm  pose  (15  poses)  and  the 
location  of  the  tracked  tool  wrt  the  FTRAC  fiducial  (5  per  pose).  For  any  known  amount  of  artificially  added 
mis-calibration,  the  series  of  frame  transformations  was  used  to  project  the  pointer  tip  on  each  image  ( pproj )  and 
then  compared  to  its  corresponding  pseg.  Thus  to  experimentally  study  the  sensitivity  of  a  VF  system  to  C-arm 
mis-calibration,  \\pproj  —  PsegW  *s  plotted  as  a  function  of  mis-calibration  in  various  graphs  as  shown  in  Figure  9. 

It  can  be  seen  that  the  overall  average  accuracy  of  out  VF  system  without  any  mis-calibration  is  about 
2.05  mm  (STD  1  mm)  in  the  2D  X-ray  images.  This  indicates  a  3D  targeting  error  of  about  1  mm  at  our  mag¬ 
nification  (0.4),  which  is  reasonable  owing  to  the  accuracies  of  (a)  the  polaris  (0.4  mm);  (b)  FTRAC  (0.6  mm); 
&  (c)  pointer  calibration  (0.2  mm). 

Figure  9  (a)  -  (c)  plot  the  mean/std  of  the  VF  system  as  a  function  of  uniformly  distributed  mis-calibration. 
Each  graph  is  computed  from  a  total  of  40,500  random  points  (75  images,  12  intervals,  45  random  runs  per 
image  per  interval).  Note  that  the  three  plots  remain  extremely  stable  &  indicate  a  similar  pattern.  The  mis- 
calibration  in  focal  length  is  slightly  more  critical  than  the  ones  in  the  origin.  Furthermore,  the  graphs  indicate 
that  a  mis-calibration,  as  high  as  even  100  mm,  does  not  result  in  any  significant  loss  in  accuracy.  To  decouple 
the  mis-calibration  errors  with  any  intrinsic  error  in  our  VF  system,  in  Figure  9  (d)  we  plot  the  distance  of 
the  mis-calibrated  projection  to  the  projection  obtained  by  using  perfect  calibration.  Thus  offsetting  for  the 
inherent  errors  in  the  system,  it  helps  us  understand  the  variation  in  the  tool  projection,  purely  as  a  function 
of  C-arm  mis-calibration.  It  can  be  observed  that  on  an  average,  the  tool  tip  projection  moves  by  less  than  a 
pixel  (0.44  mm)  for  mis-calibration  as  high  as  even  100  mm.  Figure  9  (e)  plots  the  error  in  the  VF  system  as 
a  function  of  the  distance  between  the  tool  tip  and  the  FTRAC  fiducial  used  for  the  registration.  As  indicated 
by  Equation  20,  there  is  a  slight  increase  in  the  error  of  the  mis-calibrated  VF  system  as  the  distance  between 
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Figure  9.  The  error  in  a  virtual  fluoroscopy  system  using  a  total  of  40,500  samples  plotted  as  a  function  of  mis-calibration 
in  (a)  C-arm  origin  ;  (b)  C-arm  focal  length  ;  (c)  origin  &;  focal  length  together.  It  can  be  noticed  that  VF  guidance  is 
practically  independent  of  mis-calibration.  (d)  The  additional  error  in  the  VF  system  is  plotted  as  a  comparison  to  the 
projection  with  the  perfect  calibration  at  each  image.  It  can  be  noticed  that  the  additional  error  due  to  mis-calibration  is 
not  significant,  indicating  that  the  primary  sources  of  error  are  from  other  sources  like  tool  tracking  or  tracker-fluoroscope 
registration,  (e)  The  error  in  the  VF  system  plotted  as  a  function  of  the  distance  between  the  anatomy  (tool-tip)  from 
the  registration  fiducial.  Typical  patient  thickness  of  300  mm  does  not  add  any  significant  error  to  the  system. 


the  tool-tip  and  the  fiducial  goes  over  300  mm.  Nevertheless,  this  increase  is  only  marginal  (less  than  a  pixel), 
indicating  a  significant  amount  of  robustness  to  the  tool-tip  position.  Both  graphs  in  Figure  9  (d)-(e)  use  the 
same  data  as  used  in  (c).  Thus,  it  can  be  concluded  that  mis-calibration  as  high  as  50  mm  does  not  add  any 
significant  amount  of  additional  error  in  a  VF  system.  Moreover,  this  additional  error  can  tolerate  a  350  mm 
distance  from  the  tracked-tool  to  the  registration  fiducial,  a  number  representing  typical  patient  sizes. 

Distortion  Correction:  With  the  distortion  correction  phantom  attached  to  the  intensifier,  86  C-arm  images 
in  a  25°  cone  around  the  AP-axis  were  acquired  at  a  5°  separation.  The  average  distortion  in  the  images  was 
found  to  be  3.31  mm  with  a  standard  deviation  of  1.41  mm.  Figure  10  (a)  plots  the  average  variation  in 
the  image  distortion  from  one  pose  to  another,  as  a  function  of  angular  image  separation.  Based  on  the  two 
methods  suggested  in  Section  2.3,  Figure  10  (b),  (c)  plot  the  performance  of  each  algorithm  as  a  function  of  the 
repeatability  of  C-arm  pose.  The  center-image  distortion  correction  function  reduces  the  average  distortion  in 
the  image  to  an  average  value  of  0.63  mm  (STD  0.33  mm),  using  just  a  single  center  image  in  a  25°  cone.  The 
mean-image  distortion  correction  function  performs  better,  reducing  the  residual  distortion  error  to  an  average 
of  0.48  mm  (STD  0.3  mm).  Though  the  center-image  correction  method  is  less  accurate  than  the  mean-image 
correction  method,  it  has  the  convenience  of  using  just  a  single  pre-operative  image  near  the  intended  region 
of  C-arm  use.  Thus  the  results  indicate  that  using  a  single  distortion  correction  image  can  reduce  the  average 
distortion  error  in  the  image  from  3.31  mm  to  0.63  mm  when  the  intra-operative  C-arm  pose  can  be  repeatability 
positioned  with  25°  accuracy,  which  is  sufficiently  acceptable  for  many  applications.  The  use  of  a  mean-image  by 
utilizing  multiple  images  close  to  the  intended  region  of  use,  can  further  decrease  this  error  to  0.48  mm,  though 


at  the  expense  of  collecting  multiple  images  in  the  region.  Note  that  these  numbers  for  residual  distortion  errors 
are  close  to  the  resolution  of  the  X-ray  images  (0.44  mm). 


Distortion  difference  vs.  Image  separation  Distortion  Error  from  center  image  Distortion  Error  from  mean  distortion  map 


Figure  10.  (a)  The  variation  in  C-arm  distortion  from  one  pose  to  the  other  as  a  function  of  angular  separation.  We 
notice  that  distortion  is  continuous  and  does  not  change  radically  near  a  certain  neighborhood,  (b)  The  performance  of 
using  the  geographic  center-image  for  distortion  correction  of  all  images  in  a  region.  The  residual  error  in  distortion  is 
shown  as  a  function  of  pose  repeatability  (cone  half-angle).  Only  a  single  image  is  needed  for  distortion  correction  in  the 
whole  neighborhood  of  interest,  (c)  The  performance  (residual  error)  of  using  a  computed  mean-distortion  map  for  the 
correction  of  all  images  inside  a  region.  Multiple  pre-operative  distortion  correction  images  in  the  region  will  be  needed 
to  compute  an  estimate  of  the  mean  distortion. 


4.  CONCLUSION 

We  modelled  the  the  effects  of  mis-calibration  on  as  an  affine  transform,  and  proved  its  validity  experimentally  on 
a  variety  of  common  fluoroscopy  procedures  involving  C-arm  tracking,  3D  reconstruction  or  virtual  fluoroscopy. 
We  have  derived  bounds  on  the  amount  of  scaling,  translation  and  rotation  error.  For  pose  dependant  calibra¬ 
tion,  we  proved  that  using  the  mean  calibration  minimizes  the  reconstruction  variance.  Phantom  experiments 
with  a  radiographic  fiducial  indicate  that  C-arm  tracking  is  insensitive  to  mis-calibrations.  We  also  showed  that 
mis-calibration  up  to  50  mm  adds  no  additional  error  in  3D  reconstruction  of  small  objects,  beyond  which  the 
reconstructed  objects  begin  to  drift  wrt  the  fiducial,  while  still  retaining  the  shape.  For  the  case  when  external 
trackers  are  used  in  conjunction  with  C-arms,  we  showed  that  mis-calibrations  of  up  to  50  mm  are  still  acceptable 
if  the  tool-projection  is  kept  in  the  center  of  the  image  or  if  the  tracker-fluoroscope  registration  phantom  is  kept 
close  to  the  patient  anatomy.  Experimentally,  we  showed  that  virtual  fluoroscopy  can  tolerate  mis-calibrations 
as  high  as  100  mm.  Furthermore,  to  address  the  problem  of  pose  dependant  distortion  correction,  we  propose 
the  use  of  a  mean-image,  reducing  the  distortion  close  to  1  pixel  error.  In  conclusion,  a  significant  family  of 
quantitative  fluoroscopy  applications  involving  localization  of  small  markers  can  function  without  cumbersome 
on-line  calibration.  A  constant  loose  calibration  might  suffice. 
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ABSTRACT 

For  quantitative  C-arm  fluoroscopy,  we  had  earlier  proposed  a  unified  mathematical  framework  to  tackle  the 
issues  of  pose  estimation,  correspondence  and  reconstruction,  without  the  use  of  external  trackers.  The  method 
used  randomly  distributed  unknown  points  in  the  imaging  volume,  either  naturally  present  or  induced  by  placing 
beads  on  the  patient.  These  points  were  then  inputted  to  an  algorithm  that  computed  the  3D  reconstruction.  The 
algorithm  had  an  8°  region  of  convergence,  which  in  general  could  be  considered  sufficient  for  most  applications. 
Here,  we  extend  the  earlier  algorithm  to  make  it  more  robust  and  clinically  acceptable.  We  propose  the  use  of 
a  circle/ellipse,  naturally  found  in  many  images.  We  show  that  the  projection  of  elliptic  curves  constrain  5  out 
of  the  6  degrees  of  freedom  of  the  C-arm  pose.  To  completely  recover  the  true  C-arm  pose,  we  use  constraints 
in  the  form  of  point  correspondences  between  the  images.  We  provide  an  algorithm  to  easily  obtain  a  virtual 
correspondence  across  all  the  images  and  show  that  two  correspondences  can  recover  the  true  pose  95%  of  the 
time  when  the  seeds  employed  are  separated  by  a  distance  of  40  mm.  or  greater.  Phantom  experiments  across 
three  images  indicate  a  pose  estimation  accuracy  of  1.7°  using  an  ellipse  and  two  sufficiently  separated  point 
correspondences.  Average  execution  time  in  this  case  is  130  seconds.  The  method  appears  to  be  sufficiently 
accurate  for  clinical  applications  and  does  not  require  any  significant  modification  of  clinical  protocol. 

Keywords:  Tracking,  Localization,  Registration,  X-ray  reconstruction,  C-arm  pose  estimation 

1.  INTRODUCTION 

With  an  approximate  annual  incidence  of  220,000  new  cases  and  33,000  deaths  prostate  cancer  continues  to  be 
the  most  common  cancer  in  men  in  the  United  States.1  For  several  decades,  the  definitive  treatment  for  low  risk 
prostate  cancer  was  radical  prostatectomy  or  external  beam  radiation  therapy,2  but  low  dose  rate  permanent 
seed  brachytherapy  (shortly  brachytherapy  thereafter  in  this  document)  today  can  achieve  virtually  equivalent 
outcomes.3,4  The  success  of  brachytherapy  chiefly  depends  on  our  ability  to  tailor  the  therapeutic  dose  to  the  pa¬ 
tients  individual  anatomy.  In  contemporary  practice,  however,  implant  planning  is  based  on  idealistic  preplanned 
seed  patterns  that,  as  15  years  of  clinical  practice  has  clearly  demonstrated,  are  not  achievable  in  the  actual 
human  body.  According  to  a  comprehensive  review  by  the  American  Brachytherapy  Society,5  the  preplanned 
technique  used  for  permanent  prostate  brachytherapy  has  limitations  that  may  be  overcome  by  intraoperative 
planning.  At  the  same  time,  continues  the  re-port,  the  major  current  limitation  of  intraoperative  planning  is 
the  inability  to  localize  the  seeds  in  relation  to  the  prostate.  There  are  excellent  algorithmic  and  computational 
tools  available  today  to  optimize  a  brachytherapy  treatment  plan  intraoperatively,  thereby  allowing  for  an  im¬ 
proved  dose  coverage.  These  methods,  however,  critically  require  that  the  exact  three-dimensional  locations  of 
the  implanted  seeds  are  precisely  known  with  respect  to  the  patients  anatomy. 

C-arm  fluoroscopy  is  attractive  as  a  means  to  generate  this  precise  knowledge.  It  is  the  most  widely  used 
intraoperative  imaging  modality  in  general  surgery.  However,  it  presently  lacks  the  ability  for  robust  and  easy 
quantitative  guidance. 
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Quantitative  fluoroscopy-guided  surgery  needs  to  solve  four  major  problems:  1  C-arm  image  distortion;  2  the 
calibration  of  imaging  parameters;  3  pose  recovery  or  tracking;  and  4  registration  to  imaging  modalities.  The  first 
two  are  well- studied  problems  in  the  literature.6-8  In  particular,  it  has  been  suggested  that  precise  calibration 
of  C-arm  parameters  do  not  significantly  improve  pose  reconstruction.9  Through  modeling  the  effects  of  mis- 
calibration  as  an  affine  transform,  it  is  demonstrated  that  C-arm  tracking  is  insensitive  to  mis-calibrations: 
mis-calibration  of  up  to  50  mm  adds  no  additional  error  in  3D  reconstruction  of  small  objects.  A  constant 
calibration  can  therefore  be  assumed  for  all  images.  On  the  other  hand,  pose  recovery  on  unencoded  C-arm 
machines  is  a  major  technical  problem  that  presently  does  not  have  a  clinically  practical  solution  in  many  areas 
of  application. 

In  current  commercial  C-arm  fluoroscopy  surgical  navigation  systems,  pose  recovery  is  performed  through 
localizing  the  x-ray  detector  in  room  coordinates  by  some  auxiliary  optical  tracker10,11  or  electromagnetic  EM 
tracker.12  Unfortunately,  auxiliary  trackers  sometimes  become  impractical  for  various  reasons.  They  are  expen¬ 
sive  and  add  to  the  complexity  of  the  operating  room  since  they  require  an  additional  calibration  step.  Optical 
trackers  require  a  line  of  sight,  which  becomes  cumbersome  in  a  clinical  setting  and  requires  an  alteration  in 
the  standard  work  flow.  The  EM  trackers  can  successfully  overcome  this  issue,  but  become  susceptible  to  field 
distortion  from  metal  objects  like  surgical  tools  or  the  C  arm  itself,  compromising  on  accuracy. 

In  a  recent  publication,13  the  authors  delineate  the  above  problems  and  also  say  that  using  optical  trackers 
reduces  the  useful  imaging  volume  of  the  fluoroscope  and  potentially  compromises  the  achievable  accuracies. 
Despite  using  an  optical  tracker  to  track  their  surgical  tools,  they  explicitly  choose  to  not  track  the  C  arm  using 
the  tracker  but  instead  use  a  radio-opaque  fiducial.  Their  system  has  been  fairly  successful  for  various  surgeries 
and  has  been  in  clinical  use  for  the  last  four  years.  To  make  the  fiducials  feasible,  recent  publications  have 
reported  smaller  fiducials  by  compact  bead  placement.  In  one  study  an  artificially  induced  fiducial,  consisting 
of  several  lines,  seeds  and  ellipses,  was  demonstrated  to  be  capable  of  uniquely  determining  correct  C-arm  pose 
with  mean  accuracy  of  .33°  in  rotation  and  .53  mm  in  translation.14 

Solutions  that  require  only  objects  naturally  present  in  the  image  frame  in  the  surgery  room  are  more 
attractive  still.  In  many  applications,  screw/needle  ends,  implanted  surgical  markers,  special  anatomy  points 
etc.  are  naturally  present  in  the  images.  By  enforcing  the  ”  consistency”  of  these  feature  points  across  the  images, 
one  can  potentially  solve  for  all  unknown  parameters  of  calibration,  pose  recovery,  matching,  and  reconstruction 
in  one  massive  high-dimensional  non-linear  optimization.  In  particular,  in  the  case  of  brachytherapy  there  are 
many  seeds  present  that  may  be  usable  for  optimization. 

Accurate  three  dimensional  reconstruction  of  the  brachytherapy  seeds  requires  that  individual  seeds  be 
matched  between  multiple  images.  Seed  matching  can  be  formalized  as  a  combinatorial  optimization  prob¬ 
lem,  and  demonstrates  that  the  Hungarian  algorithm  can  be  applied  to  match  seeds  between  images.  It  is  proved 
that,  using  the  Hungarian  algorithm,  two  images  are  insufficient  to  uniquely  match  seeds.  Three  images  are 
necessary,  and  in  practice  more  may  be  needed  for  the  identification  of  seeds  that  are  not  visible  in  all  images.15 

Considerable  work  has  been  done  in  studying  the  usefulness  of  conic  images  in  image  reconstruction,  both  in 
segmenting  conic  curves  in  images  and  in  determining  pose  based  on  successful  segmentations.16-24  It  has  been 
demonstrated  that  images  of  conics  can  be  used  in  pose  determination.  In  particular,  it  is  possible  to  determine, 
to  a  single  rotational  degree  of  freedom,  the  pose  from  which  an  image  was  generated  based  on  a  single  image  of 
a  circle  within  the  image  frame.25 

In  this  paper,  we  propose  a  solution  to  accurately  estimate  the  pose  by  placing  a  wire  in  the  shape  of  an  elliptic 
curve  in  the  imaging  frame.  We  demonstrate  that  images  of  this  ellipse  constrain  5  of  the  6  pose  reconstruction 
degrees  of  freedom.  Previously  designed  algorithms  that  determine  a  correspondence  between  brachytherapy 
seeds  can  be  used  to  provide  seed  correspondences.  These  correspondences  will  be  combined  with  a  non-linear 
optimization  to  determine  C-arm  pose. 


2.  METHODS 

2.1.  Imaging  Setup 

A  C-arm  fluoroscope  was  used  to  generate  three  images  of  an  artificial  brachytherapy  seed  cloud  to  be  used 
to  test  our  ability  to  reconstruct  C-arm  pose.  A  fiducial  device  that  includes  two  wire  ellipses  is  placed  in  the 


Figure  1.  Fiducial  device  containing  wire  ellipses,  and  seed  cloud  phantom  to  be  used  in  assessing  algorithm.  Within 
the  phantom  cloud  block  are  precisely  positioned  sample  seeds. 


Figure  2.  Three  sample  C-arm  Images  of  a  phantom  brachytherapy  seed  cloud  and  fiducial  with  ellipse 


imaging  frame,  and  is  used  to  create  an  image  of  an  ellipse  in  each  C-arm  image.  The  design  for  the  fiducial  was 
refined  by  generating  multiple  rough  prototypes  from  ABS  acrylonitrile  butadiene  styrene  using  an  FDM  fused 
deposition  modeling  rapid  prototype  machine.  The  final  fiducial  design  was  then  fabricated  from  a  acetol  rod 
using  a  four- axis  CNC  computer  numerical  control  mill.  The  cylinder  was  press  fit  into  a  custom  acetol  mount 
that  provided  three  mutually  orthogonal  sets  of  mounting  holes  for  attaching  to  an  accurate  rotary  table  for 
validation.  The  dimensions  of  the  ellipse  were  precisely  controlled  in  its  manufacture. 

The  brachytherapy  seed  cloud  was  modeled  by  a  highly  precise  phantom.  Real  x-ray  images  were  taken 
of  the  phantom  using  a  fluoroscope  Philips  BV  3000.  The  system-supplied  parameters  were  read  from  the 
DICOM  header,  otherwise  the  fluoroscope  was  not  explicitly  calibrated.  Moreover,  the  images  were  not  distortion 
corrected  distortion  2  mm.  The  fiducial  was  mounted  on  a  highly  accurate  0.002°  resolution  rotational  turntable 
30000  Heavy  Load  Worm  Gear  Drive  from  Parker  Automation,  Irwin,  Pennsylvania.  The  fluoroscope  remains 
stationary  while  the  fiducial  moves  in  a  known  path,  providing  ground  truth.  The  fiducial  with  seed  cloud 
phantom  attached  is  shown  in  Figure  1  A  sample  set  of  three  images  produced  by  this  set-up  is  presented  in 
Figure  2. 

A  fiducial  mount  was  designed  such  that  it  produced  zero  translation  and  a  known  rotation  when  the  turntable 
was  rotated.  The  design  supported  the  two  independent  rotation  axes  typical  to  C  arms.  Thus,  given  any  starting 
pose,  the  relative  motion  between  the  current  pose  and  the  starting  pose  is  known  precisely  from  the  turntable 
reading.  The  relative  motion  is  also  calculated  using  our  algorithm  from  computed  current  pose  and  the  computed 
starting  pose.  The  difference  between  the  computed  relative  motion  and  known  relative  motion  is  the  error. 

2.2.  Reconstruction  Algorithm  Design 

An  algorithm  was  designed  and  executed  on  the  C-arm  images  to  reconstruct  pose.  Previously  designed  segmen¬ 
tation  algorithms  were  used  to  determine  the  location  of  the  ellipse  and  the  brachytherapy  seeds  in  the  images. 
These  algorithms  produced  ellipses  equations  of  the  form: 


Au 2  +  Buv  +  Cu 2  +  Du  +  FT  +  G  =  0 


(1) 


in  each  image,  in  which  u  and  v  are  respectively  the  coordinates  along  the  x  and  y  axis  of  the  image,  measured 
in  pixels,  and  A,  B ,  C,  D,  F,  and  G  are  constants  determined  by  the  algorithm.  This  equation  and  the  assumed 
C-arm  calibration  can  be  used  to  determine  the  equation  of  the  three  dimensional  cone  which  has  as  vertex  the 
focal  point  of  the  X-ray  source.  This  equation  will  be 


A0x2  +  B0xy  +  C0x 2  +  ^ xz  +  T yz  +  ^pzz2  =  0.  (2) 
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where  x,  y,  z  are  coordinates  in  a  system  that  has  the  focal  point  of  the  x-ray  device  as  origin  with  £  axis 
perpendicular  to  the  image  plane  and  x  and  y  axis  parallel  to  the  imaging  plane,  and  /  is  the  focal  length  of  the 
imaging  scenario.  This  equation  can  be  represented  in  matrix  form  as 
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Given  the  dimensions  of  the  ellipse,  our  task  was  to  find  the  transformation  between  a  coordinate  system 
centered  on  the  physical  ellipse  with  x  axis  along  its  major  axis  and  y  axis  along  its  minor  axis  and  the  system 
with  origin  at  the  imaging  focal  point  and  with  x  axis  along  the  image  x  axis  of  the  image,  y  axis  along  the 
image  y  axis  of  the  image,  and  z  axis  perpendicular  to  the  image  plane. 

We  broke  the  process  into  several  steps  using  a  method  similar  to  that  described  for  the  case  of  a  circle  by 
Costa  and  Shapiro.25  The  first  step  involved  finding  the  rotation  between  the  system  of  coordinates  defined  by 
the  image  plane,  and  the  system  with  x  and  y  axes  along  the  major  and  minor  axis  of  the  cone  and  2  axis  along 
the  cone’s  center.  This  coordinate  system  is  of  particular  interest  because  in  it  the  cone  equation  is  of  the  form: 
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Having  the  equation  in  this  form  makes  it  easier  to  work  with  in  the  next  step.  The  rotation  matrix  that 
achieves  it  is  just  the  eigenvectors  of  the  original  cone  matrix  Mq,  since  multiplying  a  matrix  by  one  of  its 
eigenvectors  and  right  multiplying  by  the  transpose  of  that  same  matrix  produced  a  diagonalized  version  of  the 
initial  matrix. 

We  next  sought  to  find  the  transformation  between  this  coordinate  system  with  2  axis  perpendicular  to 
the  cone  base  and  one  with  2  axis  perpendicular  to  the  surface  of  the  physical  wire  ellipse.  Through  applying 
constraints  imposed  by  the  ellipse,  we  sought  to  find  the  set  of  transformations  that  would  achieve  this  for  a 
given  image.  We  parameterized  the  problem  by  defining  two  values,  and  0  which  together  can  represent  a 
rotation  from  the  frame  with  axis  aligned  with  the  cone  base  to  any  other  axis  system.  The  coordinate  system 
is  first  rotated  by  0  about  the  2  axis  of  the  original  coordinate  system,  then  by  0  around  the  new  y  axis.  This 
relationship  is  demonstrated  in  Figure  3  For  any  image,  there  are  infinitely  many  possible  <f>  —  0  combinations, 
corresponding  to  infinitely  many  possible  physical  ellipses  that  could  generate  a  given  ellipse  image.  These 
rotations  can  be  expressed  by  the  following  matrices: 
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Figure  3.  Graphical  representation  of  the  angles  phi  and  theta,  used  in  pose  reconstruction 
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Our  goal  therefore  is  to  solve  for  </>  and  0  such  that: 


(6) 


XR^ReMPRlR^XT  =  0 


(7) 


Where  Mp  is  the  equation  for  the  cone  in  the  coordinate  system  with  z  axis  perpendicular  to  the  actual 
ellipse.  We  have  several  constraints  related  to  this  cone  that  we  can  put  in  place  to  determine  the  set  of  <p  —  0 
solutions.  We  know  that  any  plane  parallel  to  the  x-y  plane  will  slice  the  cone  in  a  section  with  eccentricity  equal 
to  the  eccentricity  of  the  physical  ellipse,  and  this  eccentricity  is  a  known  constant.  So,  if  we  set  Z  constant  we 
will  get  the  equation  of  an  ellipse  with  eccentricity.  Doing  this  gives  us: 


Ap  *  x2  +  Bp  *  xy  +  Cp  *y  =  K 


(8) 


AP  =  Go  *  sin2{0)  +  Go  *  cos2 (6) sin2 ty)  +  A0  *  cos2 (0) cos2 ty) 


(9) 


Bp  =  2{— Aq  *  cos (0) cos ((f))  *  sinty)  +  Go  *  cos (0) sin (0) cos (0)} 


(10) 


Cp  =  Aq  *  sin2  (0)  +  Go  *  cos2((/)) 

The  ellipse  generated  in  this  slice  will  be  described  by  the  equation: 
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(12) 


With  i/j  as  some  rotation  angle,  and  a  and  b  as  the  major  and  minor  axes  of  the  real  ellipse.  K  will  be  some 
constant  dependent  on  the  z  value  at  which  the  slice  was  taken  Expanding  this  equation  gives  us: 


cos2ty)  sin2ty)\  2  /  1  1\  f  sin2  ty)  cos2ty)\  2 


(13) 


So,  we  know  that  for  some  angle  of  rotation  -0: 


(  cos2( xp)  sin2(!p)\ 

af={-^+—^) 

BP  =  sin(2ip)  (A  -  )-) 

(  sin2( ip)  cos2{^ip)\ 

Cp=(^^  +  ^^J 

We  can  combine  these  statements  to  place  a  constraint  on  the  allowed  values  of  these  constants  relative  to 
each  other:  trigonometric  identities  allow  us  to  say  that: 


A-+c-  =  K(h  +  h) 

(Ap  —  Cp)2  +  Bq  =  K2 

AP  +  CP  H K  +  A) 

These  can  be  combined  to  say  that: 

(ap  -  Cp)2  +  gy  (^-i)2 
( Ap  +  Cp )2  (^-  +  1)2 

Substituting  the  values  for  Ap,  Bp,  and  Cp  previously  presented  in  equations  9,  10,  11  allows  us  to  obtain 
an  equation  relating  6  to  <p>.  This  equation  is  a  quadratic  equation  in  cos2((/))  of  the  form: 

0  =  Ki{9)  *  COS4(0)  +  K2(9)  *  cos2((t>)  +  K3(9)  (14) 

With  the  Kn(6)  functions  of  0  defined  as  follows: 

Kx=-  4  *  (-Cp  +  Ap)2  *  ^  *  (1  -  2  *  cos(0)2  +  cos(ef)/( ^  +  l)2  (15) 
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Since  14  is  a  quadratic  equation  in  cos2(0),  for  any  given  value  of  0  there  are  4  values  of  9  that  satisfy  the 
constraint  that  it  imposes.  After  choosing  a  —  0  combination  to  consider,  we  apply  these  transformations  to 
our  coordinate  system.  This  gives  us  a  cone  equation  of  the  form: 

XR4)RQMPRT^RTQXT  =  XCeXt  =  0  (18) 

Where  CE  is  the  cone  in  a  coordinate  system  with  2  axis  perpendicular  to  the  face  of  the  real  ellipse.  We 
next  apply  a  rotation  by  the  angle  0  to  bring  the  x  and  y  axes  of  our  coordinate  system  parallel  to  the  major 
and  minor  axes  of  the  real  ellipse.  To  solve  for  0,  we  solve  for  the  values  that  sets  BE  to  0.  This  angle  is  given 
by: 


'0  =  1/2*  atan(2  *  BEj (C E  —  AE')') 


(19) 


Transforming  the  coordinate  system  by  a  rotation  around  the  z  axis  of  the  angle  0,  which  is  described  by  the 
transformation  matrix: 
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The  coordinate  system  created  by  the  combination  of  these  rotations  now  has  axes  parallel  to  the  axes  of  the 
real  ellipse,  but  is  centered  at  the  imaging  focal  point.  The  final  step  in  our  algorithm  is  therefore  a  translation 
to  bring  the  center  of  the  coordinate  system  to  the  center  of  the  ellipse.  This  translation  is  computed  by  taking 
a  slice  of  the  cone  a  distance  /  away  from  the  focal  points,  computing  its  major  axis  at  this  point,  and  using 
this  computation  to  determine  where  on  the  cone  a  slice  would  have  major  axis  equal  to  the  major  axis  of  the 
actual  real  ellipse.  The  coordinate  system  is  then  translated  by  this  vector. 

Combining  these  transformations  allows  us  to  create  an  algorithm  to  determine  the  set  of  possible  poses  that 
could  have  created  a  given  image.  However,  as  infinitely  many  possible  (j)  —  9  combinations  exist,  the  ellipse  image 
can  not  be  used  to  uniquely  determine  the  pose.  It  constrains  5  of  the  6  degrees  of  freedom  of  pose,  allowing  a 
single  degree  of  freedom  that  must  be  removed  through  further  constraints. 

To  constrain  this  set  to  a  single  solution,  additional  information  is  necessary.  Correspondences  between 
points  in  the  brachytherapy  seed  cloud  can  be  used  to  create  the  necessary  additional  constraints.  We  designed 
an  algorithm  that  could  take  as  input  a  proposed  set  of  poses  and  a  set  of  corresponding  point  images,  and 
determine  the  actual  three-dimensional  location  of  the  image  points,  and  the  error  in  point  segmentation  that 
this  set  of  poses  would  imply.  By  optimizing  for  minimal  projection  error,  we  were  able  to  determine  possible 
poses  that  might  produce  the  images  observed. 


3.  RESULTS 


The  algorithm  was  tested  using  C-arm  images  of  a  phantom  brachytherapy  seed  cloud  near  a  metal  ellipse. 
Initially,  trials  were  run  with  only  one  additional  point  correspondence.  This  seemed  like  an  attractive  possibility 
since  this  correspondence  could  come  from  the  mean  location  of  all  segmented  seeds  and  might  allow  for  the  seed 
correspondence  step  to  be  bypassed.  With  only  one  additional  point  correspondence,  the  algorithm  succeeded 
in  only  9  of  20  cases,  indicating  that  a  single  correspondence  is  insufficient  to  satisfactorily  resolve  pose.  Sets 
of  2  and  3  seeds  were  selected  at  random,  and  were  used  in  the  optimization  function  to  attempt  to  reconstruct 
C-arm  pose.  Average  distance  between  the  points  used  in  optimization  was  measured  and  recorded.  A  sample 
reconstruction  is  presented  in  Figure  4. 

With  the  addition  of  other  point  correspondence,  the  success  of  the  algorithm  was  shown  to  depend  con¬ 
siderably  on  the  distance  between  the  two  points  being  used.  The  algorithm  was  defined  to  have  succeeded  if 


Figure  4.  Sample  representation  of  the  pose  reconstruction  from  the  images  presented  in  Figure  2.  On  left  is  the  ground 
truth,  while  on  right  is  the  reconstructed  pose.  Each  of  the  cone  vertices  represents  an  x-ray  source,  while  the  cone 
bases  represent  ellipse  images.  The  three  cones  intersect  at  the  location  of  the  real  ellipse.  The  points  are  points  used  to 
generate  additional  constraints;  their  images  lie  in  plane  with  the  ellipses  to  which  they  correspond 
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Figure  5.  Observed  rate  of  failure.  On  the  left  are  the  results  from  optimizing  using  2  brachytherapy  seeds,  and  on  the 
right  are  results  produced  using  3  brachytherapy  seeds.  X  axis  is  the  average  spacing  between  seeds  used.spacings 
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Figure  6.  Observed  Angular  Error  in  non- failure  cases  for  for  2  and  3  brachytherapy  seeds,  using  sets  of  seeds  with 
various  average  spacings 


the  average  error  in  reconstructed  pose  as  compared  with  actual  pose  was  less  than  7°  The  rate  at  which  the 
algorithm  was  able  to  successfully  reconstruct  pose  is  presented  as  a  function  of  bead  spacing  in  Figure  5.  As 
this  figure  demonstrates,  for  most  cases  with  2  points,  or  for  cases  with  3  points  and  wide  spacing,  failure  rate 
becomes  quite  low.  For  those  cases  in  which  pose  was  reconstructed,  the  average  angular  error  still  did  display 
some  dependence  on  seed  spacing  and  number  of  seeds  used  in  reconstruction.  This  relationship  is  presented  in 
Figure  6. 


4.  DISCUSSION 

An  algorithm  has  been  presented  that  is  capable  of  determining  C-arm  pose  from  three  images  of  a  single  ellipse 
and  correspondences  between  additional  points.  Coupled  with  previous  work  that  enables  point  correspondence 
matching,  this  algorithm  could  be  applied  to  reconstruct  C-arm  pose  for  three  images  of  a  brachytherapy  patient, 
with  a  single  artificially  manufactured  wire  ellipse  of  known  dimension  as  the  only  addition  to  the  imaging 
situation.  This  represents  a  significant  step  toward  a  complete  pose  reconstruction  system  employing  no  external 
hardware. 

An  aim  of  further  work  might  be  to  determine  the  number  of  point  correspondences  that  should  be  theoreti¬ 
cally  required  to  exactly  determine  pose.  In  addition,  it  would  be  desirable  to  include  the  ellipse  into  the  error 
calculation  process;  further  work  could  focus  on  determining  the  effective  weightings  that  should  be  given  to  the 
points  and  to  the  ellipse  in  such  an  error  calculation. 
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ABSTRACT 

C-arm  images  suffer  from  pose  dependant  distortion,  which  needs  to  be  corrected  for  intra-operative  quantitative 
3D  surgical  guidance.  Several  distortion  correction  techniques  have  been  proposed  in  the  literature,  the  current 
state  of  art  using  a  dense  grid  pattern  rigidly  attached  to  the  detector.  These  methods  become  cumbersome 
for  intra-operative  use,  such  as  3D  reconstruction,  since  the  grid  pattern  interferes  with  patient  anatomy.  The 
primary  contribution  of  this  paper  is  a  framework  to  statistically  analyze  the  distortion  pattern  which  enables 
us  to  study  alternate  intra-operative  distortion  correction  methods.  In  particular,  we  propose  a  new  phantom 
that  uses  very  few  BBs,  and  yet  accurately  corrects  for  distortion. 

The  high  dimensional  space  of  distortion  pattern  can  be  effectively  characterized  by  principal  component  analysis 
(PC A).  The  analysis  shows  that  only  first  three  eigen  modes  are  significant  and  capture  about  99%  of  the 
variation.  Phantom  experiments  indicate  that  distortion  map  can  be  recovered  up  to  an  average  accuracy  of 
less  than  0.1  mm/pixel  with  these  three  modes.  With  this  prior  statistical  knowledge,  a  subset  of  BBs  can 
be  sufficient  to  recover  the  distortion  map  accurately.  Phantom  experiments  indicate  that  as  few  as  15  BBs 
can  recover  distortion  with  average  error  of  0.17  mm/pixel,  accuracy  sufficient  for  most  clinical  applications. 
These  BBs  can  be  arranged  on  the  periphery  of  the  C-arm  detector,  minimizing  the  interference  with  patient 
anatomy  and  hence  allowing  the  grid  to  remain  attached  to  the  detector  permanently.  The  proposed  method 
is  fast,  economical,  and  C-arm  independent,  potentially  boosting  the  clinical  viability  of  applications  such  as 
quantitative  3D  fluoroscopic  reconstruction. 

Keywords:  distortion,  intra-operative  distortion  correction,  principal  component  analysis,  eigen  modes 

1.  INTRODUCTION 

C-arm  fluoroscopy  is  the  most  commonly  used  intra-operative  imaging  modality  because  of  its  low  cost  and  ease 
of  use.  These  fluoroscopic  images  are  distorted  and  need  to  be  corrected  for  use  in  quantitative  intra-operative 
analysis.  This  distortion  is  caused  mainly  due  to  the  curved  nature  of  the  detector  resulting  in  radial  or  pin 
cushion  type  of  distortion  and  due  to  the  earth’s  magnetic  field  deflecting  the  electrons  resulting  in  an  S-curve 
distortion.  Boone  et  al  and  Gutierrez  et  al  discuss  the  actual  physical  phenomenon  that  causes  distortion. 
Moreover,  C-arm  distortion  is  pose-dependant,  i.e.,  it  varies  with  the  position  of  the  C-arm. 

Several  C-arm  distortion  correction  techniques  have  been  proposed  in  literature1  7  dating  back  to  90 ’s.  Typically, 
a  dense  grid  pattern  phantom  with  either  metal  beads,1,3  6  or  grooves2  is  attached  rigidly  to  the  image  intensifier 
and  is  imaged  at  different  C-arm  poses.  A  polynomial,  usually  a  two-dimensional  fourth  or  fifth  order  polynomial 
depending  on  the  number  of  beads  in  the  phantom,  is  used  to  model  the  distortion  by  mapping  the  bead 
coordinates  in  the  fluoroscopic  phantom  image  to  the  correct  geometrical  coordinates  of  the  phantom.  An 
interpolation  function  is  then  used  to  define  the  distortion  map  for  all  the  pixels  in  the  image.  These  techniques 
are  robust  and  are  very  effective  for  offline  distortion  correction. 

These  grid-based  methods  become  cumbersome  for  intra-operative  use  (such  as  for  3D  reconstruction),  since 
the  grid  pattern  either  significantly  corrupts  the  patient’s  image,  or  necessitates  two  X-ray  shots  at  each  pose 
(one  with  and  the  other  without  the  grid).  Fahrig  et  al  proposes  to  track  C-arm,  create  a  look  up  table  for 
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a  distortion  map  at  each  pose  and  then  interpolate  the  near-by  poses  to  determine  the  distortion  map  intra- 
operatively.  This  method  is  fast  and  overcomes  the  above  mentioned  issues,  but  requires  the  C-arm  to  be  encoded 
accurately.  The  other  option  would  be  to  remember  the  poses  for  which  the  distortion  is  measured  offline  and 
then  repeat  the  procedure  intra-operatively.2  This  method  requires  that  the  C-arm  is  tracked  and  the  pose  is 
repeatable.  Intra-operative  C-arm  tracking  could  be  slow  and  expensive. 

In  this  paper,  we  propose  a  statistical  framework  which  enables  us  to  analyze  the  distortion  maps  at  different 
poses  and  extract  any  inherent  patterns.  We  have  performed  principal  component  analysis  to  create  a  statistical 
model  for  distortion  patterns  exhibited  by  a  C-arm.  Although  C-arm  distortion  patterns  appear  to  be  random 
and  complicated,  our  analysis  shows  that  they  can  be  predicted  from  the  principal  modes  of  variation  with  a 
very  decent  accuracy.  In  particular,  we  present  a  novel  intra-operative  distortion  correction  application  with 
a  phantom  which  uses  very  few  beads  on  the  periphery  of  the  C-arm  image  intensifier  along  with  the  prior 
statistical  knowledge.  This  method  is  fast,  economical,  and  C-arm  independent,  potentially  boosting  the  clinical 
viability  of  applications  such  as  quantitative  3D  fluoroscopic  reconstruction.  The  motivation  for  this  work 
(effective  distortion  correction)  might  become  obsolete  in  case  of  flat  panel  detector  C-arms  which  would  replace 
existing  C-arms  in  future.  However,  we  believe  that  this  work  presents  a  statistical  framework  to  characterize 
the  distortion  patterns  effectively  and  to  build  applications  based  on  this  analysis. 

2.  STATISTICAL  CHARACTERIZATION  OF  C-ARM  DISTORTION 
2.1.  Method 

We  have  designed  a  grid  phantom  with  steel  beads,  approximately  300  beads,  placed  at  1  cm  distance  between 
beads  in  a  square  pattern.  Several  fluoroscopic  phantom  images  were  taken  in  various  C-arm  poses  by  rigidly 
attaching  the  phantom  to  the  image  intensifier.  A  basic  image  processing  algorithm  is  used  to  segment  the  bead 
positions  in  all  the  images.  These  segmented  bead  positions  are  referred  to  as  image  coordinates  (ud,vd)  •  The 
actual/nominal  bead  positions,  (uQ,  vQ),  are  determined  from  the  physical  geometry  of  the  grid  phantom  and  the 
distortion  is  modeled  using  a  fifth  order  Bernstein  polynomial  as  shown  in  the  following  equation. 


where 
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For  each  image,  a  bezier  matrix  B  is  constructed  using  the  nominal  coordinates  of  the  bead  positions  and 
the  above  equation  can  be  written  in  matrix  form  as  follows: 

{U0bs)mx  2  =  (B(Un0rn)')rnx(n+ 1)2  *  C\n+  l)2x2  (2) 

where  U0i,s  is  the  array  of  distorted  bead  coordinates  in  the  image;  Unorn  represents  the  nominal  coordinates 
of  the  beads;  and  m  being  number  of  grid  points  on  the  phantom  (here  m  =  300). 

Define  a  distortion  vector  as 


A  d  =  (A u,  Av)  =  (ud,  vd)  -  (uo,  v0) 

The  coefficient  matrix  C  can  be  computed  in  least  squares  sense  which  minimizes  the  error  A  d  for  each  grid 
point.  This  polynomial  effectively  models  the  distortion.  After  computing  the  polynomial,  the  distortion  vectors 
for  each  pixel  in  the  distorted  image  is  computed  using  2  and  the  image  is  sampled  in  the  resulting  distortion 
map  to  create  a  rectified  image.  In  practice,  this  technique  has  proved  to  be  robust  and  accurate  with  sub-pixel 
accuracy. 


This  procedure  is  repeated  for  all  the  images  acquired  at  different  C-arm  poses  to  create  distortion  maps. 
In  order  to  extract  patterns  from  these  distortion  maps,  we  have  performed  principal  component  analysis.  The 
data  matrix  for  this  analysis  is  created  by  vectorizing  the  image  matrix  (distortion  map)  and  by  stacking  each 
vector  into  a  matrix  A  as  shown  below: 
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where  ( Auij,Avij )  is  the  distortion  vector  of  the  ith  pixel  in  the  jth  image;  m  is  the  number  of  pixels  in  the 
image;  and  k  is  the  number  of  images. 

For  computational  issues,  we  have  ignored  the  background  pixels  and  have  only  considered  the  pixels  in  the 
active  area  of  the  C-arm  detector  while  creating  the  data  matrix.  The  mean  distortion  map  is  computed  by 
taking  average  of  all  the  distortion  maps.  This  mean  value  is  then  substracted  from  the  data  matrix  to  make 
it  a  zero  mean  data  matrix,  which  is  essential  for  principal  component  analysis.  Principal  modes  of  variation 
(eigen  modes)  are  then  extracted  from  the  covariance  matrix  of  the  input  data  matrix.8  From  this  analysis,  we 
can  express  any  distortion  map  as  a  linear  combination  of  the  principal  modes. 
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where  Mo  is  the  mean  distortion  map  and  D{  is  the  ith  principal  component.  The  orthogonal  eigen  modes 
represent  new  basis  for  the  subspace  in  which  the  C-arm  distortion  patterns  exist  as  compared  to  the  very 
high  dimensional  parameterization  using  polynomials.  For  example,  a  two  dimensional  fifth  order  Bernstein 
polynomial  has  72  coefficients. 

2.2.  Results 

We  have  performed  principal  component  analysis  on  two  sets  of  distortion  maps  extracted  from  the  grid  phantom 
images.  For  the  first  dataset,  we  have  collected  120  images  of  the  phantom,  approximately  one  every  3  degrees 
along  the  propeller  axis  of  an  OEC  9600  C-arm.  We  have  relied  on  the  gradation  marks  on  the  C-arm  to  position 
the  C-arm  and  did  not  use  any  trackers/encoders.  This  would  be  a  typical  dataset  for  3D  reconstruction  with  2D 
fluoroscopic  images.  For  the  OEC  9600  C-arm,  there  are  two  data  outputs:  1)  DICOM  output  through  floppy 
and  2)video  output  through  a  frame  grabber.  Typically,  floppy  output  is  slow  but  accurate  as  compared  to  the 
fast  but  noisy  video  signal.  The  image  size  is  459x484  pixels  and  480x720  pixels  for  floppy  and  video  outputs 
respectively  and  the  pixel  size  is  0.44  after  distortion  correction.  We  have  compared  both  the  dicom  and  the 
video  images  to  verify  that  the  video  signal  does  not  affect  the  distortion.  All  these  images  were  corrected  using  a 
fifth  order  Bernstein  polynomial.  We  then  extracted  the  distortion  maps  for  each  image  and  performed  principal 
component  analysis  on  these  distortion  maps  as  explained  above. 

This  analysis  shows  that  the  first  three  principal  modes  of  variation  are  significant  and  they  capture  about 
99%  of  the  variation  as  shown  in  Figure  1(a).  Figure  1(b)  shows  the  results  of  leave-out  validation  test.  In  this 
test,  of  the  120  training  images,  we  have  used  60  images,  roughly  every  6  degrees,  to  create  a  statistical  model 
and  then  fit  this  model  to  the  remaining  60  images.  The  error  measurement  between  the  true  distortion  map 
and  the  one  estimated  from  the  model  is  given  by 
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(a)  Percentage  variance  explained  by  eigen  modes 


(b)  Residual  error  from  leave-out  validation  test 


Figure  1.  Principal  component  analysis  on  distortion  maps,  (a)  Percentage  variation  explained  by  the  dominant  eigen 
modes.  Note  that  the  first  three  modes  capture  about  99%  of  the  variation,  (b)  Residual  error  plot  after  projecting  the 
leave-out  distortion  map  onto  the  eigen  space  as  function  of  number  of  modes.  The  average  error  is  less  than  1mm  when 
the  first  three  principal  modes  are  used.  Note  that  the  results  from  dicom  and  video  images  are  very  similar. 


The  average  raw  distortion  as  measured  from  the  grid  phantom  is  about  2.1488  mm/pixel  with  a  standard 
deviation  of  1.1469  mm/pixel  and  a  maximum  distortion  of  4.6003  mm/pixel.  With  the  leave-out  validation 
test,  we  were  able  to  achieve  an  average  distortion  correction  accuracy  of  about  0.0632  mm/pixel  and  a  standard 
deviation  of  0.0294  mm/pixel  with  a  maximum  distortion  of  0.202  mm/pixel  using  the  first  three  dominant  eigen 
modes.  The  residual  errors  for  video  images  were  close  to  the  ones  from  DICOM  images  and  hence  validating 
the  fact  that  either  output  can  be  used  for  statistical  analysis.  These  results  indicate  that  eventhough  the  C- 
arm  distortions  are  complicated  and  high-dimensional,  they  can  be  represented  by  the  eigen  modes  and  can  be 
predicted  fairly  accurately. 

The  actual  distortion  modes  from  the  principal  modes  are  shown  in  Figure  2.2.  The  first  principal  mode 
appears  similar  to  barrel  distortion,  second  mode  represents  S-curve  distortion  and  the  third  mode  looks  like  a 
spiral  distortion.  We  have  predominantly  noticed  S-curve  distortion  in  all  the  images  taken  with  that  C-arm. 
While  collecting  the  data,  we  have  noticed  that  the  images  were  rotating  in  one  direction  until  halfway  and  then 
in  opposite  direction  for  the  rest  of  the  sweep.  We  came  to  know  from  the  manufacturer  that  as  we  move  the  ”C” 
through  it’s  propeller  angles,  there  is  a  shift  of  weight  from  the  bearings  on  one  side  of  the  detector  assembly,  to 
the  bearing  set  on  the  opposite  side.  This  causes  a  sudden  jump  and  also  rotates  the  images  through  the  sweep. 
We  have  included  this  as  part  of  the  distortion  model  because  there  is  no  other  way  to  explicitly  measure  these 
inconsistencies  and  to  our  surprise  these  inconsistencies  showed  up  in  the  mode  distribution  plots  for  the  mode 
weight  parameters  (Figure  2.2).  The  Figure  3(a)  clearly  shows  the  rotational  aspect  of  the  C-arm  through  the 
propeller  axis  which  is  in  agreement  with  the  first  mode  (barrel  distortion  due  to  the  curvature  of  the  detector). 
This  distortion  is  more  or  less  assumed  to  be  constant.4,9  The  rotation  caused  by  the  gantry  can  be  clearly  seen 
in  this  plot.  The  jump  motion  can  be  attributed  to  the  jumps  in  plots  given  in  Figures  3(b)  and  3(c). 

For  the  other  experiment,  we  have  collected  200  images  distributed  more  or  less  uniformly  in  the  sphere 
spanned  by  the  ”C”  of  the  C-arm.  We  have  captured  only  video  images  this  time  as  it  was  convenient  to  copy 
images  using  the  frame  grabber.  Similar  eigen  analysis  on  this  data  says  that  the  first  four  eigen  modes  capture 
99%  of  the  variation.  The  distortion  maps  defined  by  the  principal  modes  are  shown  in  Fig  2.2.  The  first  two 
modes  look  like  barrel  distortion  along  x,  y  dimensions  respectively.  The  second  mode  shows  the  S-curve  and  the 
third  mode  is  similar  to  spiral  distortion.  From  the  leave-out  validation  test,  we  are  able  to  achieve  an  average 
accuracy  of  0.07  mm/pixel,  with  0.0553  mm/pixel  standard  deviation  and  0.3578  mm/pixel  maximum  distortion, 
comparable  to  the  results  from  the  first  dataset. 

These  two  datasets  are  collected  with  the  same  phantom  on  the  same  C-arm  in  the  same  room  but  at  different 
time  periods.  The  C-arm  is  displaced  with  in  the  room  during  this  time.  Qualitatively,  the  distortion  patterns 
from  these  two  datasets  look  similar.  In  order  to  compare  the  distortion  patterns  quantitatively,  we  measured  the 


(a)  Mean  distortion  map 
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Figure  2.  Distortion  maps  as  shown  by  the  first  three  eigen  modes  from  principal  component  analysis.  The  first  mode  is 
similar  to  barrel  distortion,  second  mode  appears  close  to  an  S-curve  distortion  and  the  third  mode  shows  spiral  distortion 

projection  errors  using  equation  (4)  by  projecting  one  dataset  on  to  the  other  subspace.  This  measure  quantifies 
how  well  one  model  can  generalize  the  data  from  the  other  model.  The  average  distortion  correction  accuracy  is 
0.2329  mm/pixel,  with  a  standard  deviation  of  0.1157  mm/pixel  and  a  maximum  distortion  correction  of  0.6589 
mm/pixel.  These  results  show  that  the  statistical  model  created  from  the  volume  data  is  capable  of  modeling 
the  distortion  patterns  extracted  from  the  full  sweep  data  acquired  at  a  different  time.  These  results  strongly 
indicate  that  the  C-arm  distortion  patterns  are  not  random  and  can  be  predicted  with  reasonable  accuracies. 
A  more  controlled  experiment  with  the  same  C-arm  at  different  locations  would  provide  more  insights  into  this 
problem. 

3.  A  NOVEL  INTRA-OPERATIVE  DISTORTION  CORRECTION  APPLICATION 

Building  upon  the  statistical  framework  presented  in  the  previous  section,  we  present  a  new  intra-operative 
distortion  correction  method  with  a  phantom  that  uses  very  few  beads.  We  propose  a  new  phantom  design  with 
beads  on  the  periphery  of  the  image  intensifier.  These  beads  can  be  arranged  in  any  pattern  (circular,  square 
etc).  The  main  idea  is  to  use  the  prior  statistical  knowledge  and  the  beads  to  drive  the  distortion  correction 
process.  This  design  suits  intra-operative  needs  very  well  because  it  has  very  minimal  interference  with  patient 
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Figure  3.  Distribution  of  the  mode  weight  parameters  for  120  images  from  a  full  sweep  along  the  propeller  axis 
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Figure  4.  Distortion  maps  as  shown  by  the  first  four  eigen  modes  from  principal  component  analysis.  The  first  two 
modes  represent  barrel  distortion  in  x,y  directions,  second  mode  shows  S-curve  type  of  distortion  and  the  third  mode  is 
spiral  distortion 


anatomy  and  hence  allows  the  grid  to  remain  attached  to  the  detector  through  out  the  surgery.  The  proposed 
method  is  fast,  economical,  and  C-arm  independent,  potentially  boosting  the  clinical  viability  of  applications 
such  as  quantitative  3D  fluoroscopic  reconstruction.  Figure  3  shows  our  simulation  results  on  the  number  of 
beads  needed  to  recover  distortion  parameters  using  the  prior  knowledge.  In  this  simulation,  we  have  used 
images  from  the  full  sweep  dataset.  Of  the  120  images,  we  have  used  60  images  to  create  the  mean  distortion 
map  and  the  distortion  modes.  The  remaining  60  images  were  used  to  do  distortion  correction.  A  band  of  pixels 
was  selected  in  each  of  these  60  images,  from  which  ”n”  number  of  pixels  were  selected  as  bead  locations.  These 
pixels  are  selected  randomly  and  for  each  experiment  we  have  performed  100  trials  for  all  the  60  images.  For 
each  iteration,  we  have  optimized  mode  weights  by  minimizing  a  similarity  measure.  The  similarity  measure 
is  determined  by  the  ”  ‘squareness”  ’  of  the  bead  locations.  By  squareness,  the  bead  locations  after  correction 
should  be  on  a  line  and  the  lines  should  be  orthogonal.  In  this  optimization  step,  there  is  no  information  about 
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Figure  5.  Residual  distortion  correction  error  as  function  of  number  of  beads  used  for  correction 

the  orientation  of  the  phantom  with  respect  to  the  image  intensifier  and  hence  the  bead  locations  were  made  to 
fit  a  specific  pattern,  for  example  a  square.  We  have  used  downhill  simplex  optimization  method  in  MATLAB 
7.0  with  the  search  space  being  defined  by  the  distortion  modes  from  PC  A.  A  linear  combination  of  the  modes 
would  give  a  distortion  map.  This  distortion  map  is  interpolated  to  determine  the  rectified  locations  for  the 
beads  using  bilinear  interpolation.  The  results  indicate  that  as  few  as  14  or  15  beads  are  required  to  recover  the 
distortion  parameters  with  an  average  accuracy  of  0.17  mm/pixel  and  a  maximum  of  2.86  mm/pixel,  accuracy 
sufficient  for  most  clinical  applications.  The  maximum  error  is  high  because  the  objective  function  is  not  able 
to  capture  the  distortions.  A  slight  modification  to  the  objective  function  incorporating  the  orientation  of  the 
phantom  might  improve  the  results  significantly. 


Figure  6.  Simulation  of  the  new  phantom  design,  (a)  shows  the  fluoroscopic  image  of  a  knee  with  the  phantom  beads 
highlighted  in  red  (b)  distortion  corrected  image  using  the  dense  grid  phantom  (c)  true  distortion  map  overlaid  on 
difference  image  ((a)-  (b)).  (d)  distortion  corrected  image  using  the  beads  from  (a)  and  statistical  modes  (e)  residual 
distortion  map  after  correction  using  new  phantom  overlaid  on  the  gray  scale  difference  image  (  (c)-(d)  ) 


4.  CONCLUSIONS  AND  FUTURE  WORK 


In  this  paper,  we  have  presented  a  statistical  framework  to  analyze  C-arm  distortion  patterns.  Principal  compo¬ 
nent  analysis  on  distortion  maps  at  different  orientations  of  the  C-arm  reveal  the  inherent  patterns.  Our  results 
from  two  datasets  indicate  that  even  though  the  distortion  patterns  are  high  dimensional  and  complicated,  they 
can  be  predicted  fairly  accurately  using  the  first  three  to  four  principal  modes.  This  analysis  can  be  easily  ex¬ 
tended  to  study  the  distortion  patterns  over  time,  and  at  different  locations  and  across  different  C-arms,  adding 
more  dimensions  to  the  problem.  Such  analysis  would  give  better  insights  in  to  modeling  distortion  patterns 
effectively.  Apparently,  this  work  is  not  applicable  to  flat  panel  detector  C-arms.  Nevertheless,  any  one  using 
non- flat-panels  have  to  go  through  the  tedious  distortion  correction  procedure.  In  such  cases,  we  consider  that 
this  statistical  approach  would  address  most  of  the  current  issues  with  the  C-arm  distortion  correction. 

We  have  described  a  novel  intra-operative  distortion  correction  method  that  uses  a  phantom  with  very  few 
beads.  Our  simulation  experiments  have  shown  that  as  few  as  15  beads  can  recover  the  distortion  parameters 
with  an  average  accuracy  of  0.17  mm/pixel.  This  new  phantom  has  very  minimal  interference  with  the  patient 
anatomy  and  hence  can  remain  attached  to  the  C-arm  during  surgery.  We  believe  that  this  method  is  fast  and 
practical  and  can  be  used  effectively  for  intra-operative  applications  such  as  3D  reconstruction.  Moving  further 
away  from  the  phantoms,  another  approach  would  be  to  do  a  phantom  less  distortion  correction  of  patient  images 
using  patient  CT  as  a  fiducial.  The  2D  views  from  3D  patient  CT  can  be  compared  to  the  actual  x-ray  images 
of  the  patient  to  optimize  the  distortion  modes,  hence  catering  to  the  needs  of  intra-operative  applications. 
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